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Abstract 

these notes wore written for an introductory course on the application of multigrid 
methods to elliptic and hyperbolic partial differential equations for engineers, physicists and 
applied mathematicians. The use of more advanced mathematical tools, such as functional 
analysis, is avoided. 1 he course is intended to be accessible to a wide audience of users 
of computational methods. We restrict ourselves to finite volume and finite difference dis- 
cretization. I he basic principles are given. Smoothing methods and Fourier smoot hing 
analysis are reviewed. The fundamental multigrid algorithm is studied. The smoothing 
and coarse grid approximation properties are discussed. Multigrid schedules and structured 
programming of multigrid algorithms are treated. Robustness and efficiency are considered. 


* Research was supported by the National Aeronautics and Space Administration under NASA Contract 
No. NAS1-UM80 while the author was in residence at the Institute for Computer Applications in Science 
and Engineering (ICASE), NASA Langley Research Center. Hampton. VA 23681-0001 . 
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1 Introduction 

Readership 

The purpose of these notes is to present, at graduate level, an introduction to the application 
of mult.gHd methods to elliptic and hyperbolic partial differential equations for engineers, 
y icis s an app ed mathematicians. The reader is assumed to be familiar with the basics 
of the analysis of partial differential equations and of numerical mathematics, but the use 
more adv a n c ed mathematical tools, such as functional analysis, is avoided. The course 
is m ended to be accessible to a wide audience of users of computational methods. We do 
’ there [ or ^ ^ve deeply into the mathematical foundations. This is done in the excellent 
monograph by Hackbusch [57], which treats many aspects of multigrid, and also contains 
many practical details The book [141] is more accessible to non-mathematicians, and pays 
more attention to applications, especially in computational fluid dynamics. 

ters OH851 1 andtr 40 ^ ^ ^ “ th ® Brandt [20 ^ the first three chap- 

ters of [85] and the short elementary introduction [27], The notes are based on parts of [1411 

where further details may be found, and other subjects are discussed, notably applications in 
computational fluid dynamics. 

Significance of multigrid methods for scientific computation 

Neediess to say, elliptic and hyperbolic partial differential equations are, by and large, at the 

computations “oft Th !, V* and giving rise to extensive 

computations. Often the problems that one would like to solve exceed the capacity of even 

he most powerful computers, or the time required is too great to allow inclusion of advanced 

mathematical models in the design process of technical apparatus, from microchips to aircraft 

making design optimization more difficult. Multigrid methods are a prime source of impor- 

an a vances in a gonthmic efficiency, finding a rapidly increasing number of users. Unlike 

o her known methods, multigrid offers the possibility of solving problems with N unknowns 

with O(N) work and storage, not just for special cases, but for large classes of problems. 

Historical development of multigrid methods 

\fll b T d ° n the multig " d bibliography in [85], illustrates the rapid growth of the 
multigrid literature, a growth which has continued unabated since 1985. 

As shown by Table 1.0.1, multigrid methods have been developed only recently. In what 
probably was the first ‘true’ multigrid publication, Fedorenko [43] formulated a multigrid al- 
gon m for the standard I five-point finite difference discretization of the Poisson equation on 
a square, proving that the work required to reach a given precison is O(N). Thfcwork was 
generalized to the central difference discretization of the general Unear elUptic partial differ- 
en la equation ( . .1) in ft _ (0, 1) x (0, 1) with variable smooth coefficients by Bachvalov [8]. 
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The theoretical work estimates were pessimistic, and the method was not Ql^who 
the time The first practical results were reported in a pioneer, ng paper by Brandt [19], who 
published another paper in 1977 [20], clearly outlining the main principles and ^ Practmal 
utility of multigrid methods, which drew wide attention and marked the beg jinnmg of rapid 
development. Ihe multigrid method was discovered independently by Hackbpc 0], who 
laid firm mathematical foundations and prov.ded rebable methods ([52], [53], [54]). A 
port by Frederickson [47] describing an efficient multigrid algorithm for the Poisson equation 
led the present author to the development of a similar method for the vor 
tion formulation of the Navier- Stokes equations, resulting m an efficient method ([•■],[ ])• 

At first there was much debate and scepticism about the true merits of multigrid methods. 
Only after sufficient initiation satisfactory results could be obtained. This led a numb 
researchers to the development of stronger and more transparent convergence proofs ([4] [93] 

1141 1511 [541 [136], [137]) (see [57] for a survey of theoretical developments). Altho g 

oft vergen e oofs' of multigrid methods are complicated, their structure has now become 
more or less s.andar.ized and trasparenr. Other outhors have tried to spead confidence m 
multigrid methods by providing efficient and reliable computer programs as much , M ■ 
ble of ‘black-box’ type, for uninitiated users. A survey will be given later. The multigrid 
guide’ of Brandt ([16], [23]) was provided to give guidelines for researchers writing eir own 

multigrid programs. 


Year 64 66 71 72 73 75 76 77 

Number 1 1 1 \ \ j 3 11 


78 79 80 81 82 83 84 

10 22 31 70 78 96 94 


85 

149 


Table 1.0.1: Years number of multigrid publications 


The following topics will not be treated here: parabolic equations eigenvalue problems and 
integral equations. For an introduction to the application of mult.gr, d methods to these 
subjects, see [56], [57] and [18]. There is relatively Uttle mater, al ,n there areas although 
multigrid can be applied profitably. For important recent advances m the field o mteg al 
equations see [25] and [130]. A recent publication on parabolic multigrid is [91]. Fin 
element methods will not be discussed, but finite volume and finite drfference discretization 
will be taken as the point of departure. Although most theoretical work has been done ,n a 
variational framework, most applications us, finite volumes or finite differences The prrna- 
ples are the same, however, and the reader should have no difficulty in applying the principles 
outlined in this book in a finite element context. 
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Mu t, grid principles are much more widely applicable than jus, ,|,e numerical solution 
of differential and integral equations. Applications in such diverse areas as control theory 
optmmatioii, pattern recognition, computational tomography and particle physics are begin ’ 

|isf ° aPP ' !ar ' F ° r * Sl ' rV '‘ y ° f lhe WKl ' 1 '“"Bine applicability of imiltigrid principles, see [ 1 7 ], 


Notation 

The notation is explained as it occurs. Latin letter like u denote 
bold version u denotes a grid function, with value Uj in grid point ,r 
approximation of u(xj). 


u 11 known 
j y intended 


functions. The 
as the discrete 


2 The basic principle of multigrid methods for partial differ- 
ential equations 


2.1 Introduction 


n his chapter the basic principle of multigrid for partial differential equations will be ex- 
pained by studying a one-dimensional model problem. Of course, one-dimensional problems 
do no require application of multigrid methods, since for the algebraic systems that result 
rom discretization direct solution is efficient, but in one dimension multigrid methods can be 
analysed by elementary methods, and their essential principle is easily demonstrated. 


Introductions to the basic principles of 
[141]. More advanced expositions are given 


imiltigrid methods are given by [20], [27], [28] and 
by [1 12], [16] and [57], Chapter 2. 


2.2 The basic principle 

One-dimensional model problem 

The following model problem will be considered 

- d 2 u/dx 2 = f(x) in Q = (0, 1 ), u(0) = du(l)/dx = 0 
A computational grid is defined by 


G - {x G Bi : x = X j = jh, >=1,2, ...,2n, h = l/2n} 

The points { 2 ;^} are called the vertices of the grid. 

Equation (2.2.1) is discretized with finite differences as 

h~ 2 ( 2u, - u 2 ) = /, 
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(2.2.3) 


2«i -«,+,) = /<■ j = 2,3,...,2n- 1 

h~ 2 (-U 2 „-l + «2n) = 2^ 2n 

zSE'^^^jrarrr^^"T=« 

method. The system (2.2.3) is denoted by 


Au = f 


(2.2.4) 


In multidimensional applications of finite difference methods, the matrix A is large and sparse 
and* the nonzero pattern has a regular structure. These circumstances favour the use of 
iterative methods for solving (2.2.4). We will present one such ‘ ^ 

iterand by a superscript m, the Gmss-Sddel Herat, on method for solving (2.2.3) 

by, assuming an initial guess u° is given, 


2 u? 


= u 


m — 1 


+ h 2 h 


u'JLi + 2 u™ = u^+1 1 + h2 fj' i = 2 ’ 3, ”’ ,2n_ 1 


(2.2.5) 


A .m I rf , 771 — 

“ u 2n-l + U 2n 




For'ease of^yX wTrepuTtb s'boundary conditions by periodic boundary condition 


u(l) = u(0) ( 2 - 2 - 6 ^ 

Then the error e m = u m - u°° is periodic and satisfies 

- e™! + 2e™ = , ef = ef +2n ( 2 - 2J ) 

As will be discussed in more detail later, such a periodic grid function can be represented by 
the following Fourier series: 


e™ = C exp(ij9 a ), B a = rra/n 

a=-n+l 


(2.2.8) 


Because of the orthogonality of {e^}, it suffies to substitute ef' 1 - in (2.2.7). 

This gives ef = c™e'> e ° with 

C = g(«a)C-' , *(*.) = /< 2 - > (2 ' 2 ' 9) 
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The function g(0 a ) is called the amplification factor. It measures the growth or decay of a 
Fourier mode of the error during an iteration. We find 

\g{0 a ) | = (5 - 4cos0 o )~ 1/2 (2.2.10) 

At first it seems that Gauss-Seidel does not converge, because 

max{| 5 (0 a )| :6 a = na/n, a = -n + 1, -n + 2, = |</(0)| = 1 (2.2.11) 

However, with periodic boundary conditions the solution of (2.2.1) is determined up to a 
constant only, so that there is no need to require that the Fourier mode a = 0 decays during 
iteration. Equation (2.2.11), therefore, is not a correct measure of convergence, but the 
following quantity is: 

max{| 5 r( 0 a )| : 0 a = jra/n, a = -n + 1, -n + 2, n, a ± 0} = | 5 (^)| 

= {1 - 20? + 0(0?)} -1 / 2 = 1 - 4ir 2 h 2 + 0(h 4 ) . (2.2.12) 

It follows that the rate of convergence deteriorates as h j 0. Apart from special cases, 

in the context of elliptic equations this is found to be true of all socalled basic iterative 

methods (more on these later; well known examples are the Jacobi, Gauss-Seidel and successive 
over-relaxation methods) by which a grid function value is updated using only neighbouring 
vertices. This deterioration of rate of convergence is found to occur also with other kinds of 
boundary conditions. The purpose of multigrid is to avoid this deterioration, and to achieve 
a rate of convergence which is independent of h. 

The essential multigrid principle 

The rate of convergence of basic iterative methods can be improved with multigrid methods. 
The basic observation is that (2.2.10) shows that |</(0„)| decreases as a increases. This 
means that, although long wavellength Fourier modes (o close to 1) decay slowly (|p(0 a )| = 

1 - 0(h 2 )), short wavelength Fourier modes are reduced rapidly. The essential multigrid 
principle is to approximate the smooth (long wavelength) part of the error on coarser grids. 
The non-smooth or rough part is reduced with a small number (independent of /i) of iterations 
with a basic iterative method on the fine grid. 

Fourier smoothing analysis 

In order to be able to verify whether a basic iterative method gives a good reduction of the 
rough part of the error, the concept of roughness has to be defined precisely. 

Definition 2.2.1 The set of rough wavenumbers Q T is defined by 

0 r = {0 a = 7 ra/n, |a| > cn, a = -n + 1 , -n + 2, ..., n} (2.2.13) 
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where 0 < c < 1 is fixed constant independent of n. 

The performance of a smoothing method is measured by its smoothing factor p, defined as 
follows. 

Definition 2.2.2 The smoothing factor p is defined by 

p = max{|#(0 a )| : 0 a € 0 r } (2.2.14) 

When for a basic iterative method p < 1 is bounded away from 1 uniformly in h, we say that 
the method is a smoother. Note that p depends on the iterative method and on the problem. 
For Gauss- Seidel and the present model problem p is easily determined. Equation (2.2.10) 
shows that \g\ decreases monotonically, so that 

p = (5 — 4 cosctt ) -1 / 2 (2.2.15) 

Hence, for the present problem Gauss-Seidel is a smoother. 

It is convenient to standardize the choice of c. Only the Fourier modes that cannot be 
represented on the coarse grid need to be reduced by the basic iterative method; thus it is 
natural to let these modes constitute 0 r . We choose the coarse grid by doubling the mesh-size 
of G. The Fourier modes on this grid have wavenumbers 0 a given by (2.2.8) with n replaced 
by nf 2 (assuming for simplicity n to be even). The remaining wavenumbers are defined to 
be non-smooth, and are given by (2.2.13) with 

c = 1/2 (2.2.16) 

Equation (2.2.15) then gives the following smoothing factor for Gauss-Seidel 

p = 5~'G (2.2.17) 

This type of Fourier smoothing analysis was originally introduced by Brandt [20]. It is a 
useful and simple tool. When the boundary conditions are not periodic, its predictions are 
found to remain qualitatively correct, except in the case of singular perturbation problems, 
to be discussed later. 

With smoothly varying coefficients, experience shows that a smoother which performs well 
in the ‘frozen coefficient’ case-, will also perform well for variable coefficients. By the ‘frozen 
coefficient’ case we mean a set of constant coefficient cases, with coefficient values equal to the 
values of the variable coefficients under consideration in a sufficiently large sample of points 
in the domain. 
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Exercise 2.2.1 Determine the smoothing factor of the damped Jacobi method (defined later) 
this" sl«her dary COnditi °'’ S (2 ' 2 ' 6) ' N ° te tha ‘ With dampi "S Pieter to = i 

f E 2 Th)tuh n ;2 WIT’”!,'"' smoothin 6 °f ‘he Jacobi method applied to problem 

Note th« L L«v f ” 7 C0I, r " (0) = “ (1 > = .»■ ^ «™6 the Fourier sine series. 

mg ac or is the same as obtained with the exponential Fourier series. 

Exercise 2.2.3 Determine the smoothing factor of the Gauss-Seidel method for central 
.scretea ton of the convection-diffusion equation cdu/dx - ed*u/dx* = /. Show that for 
\c\n/£ i and c < — 1 we have no smoother. 


2.3 The two-grid algorithm 
Coarse grid approximation 

A coarse grid G is defined by doubling the mesh-size of G: 


G - {x e R :x = Xj =jh, j = 1,2, n, h = l/ n ) (2.3.1) 

The vertices of G also belong to G; thus this is called vertez-centend coarsening. The original 
grid G is called the fine grid. Let ongmai 


U :G — Bl, (J :G -» 1R 


(2.3.2) 


be the sets of fine and coarse grid functions, respectively. A prolongation operator PU^U 
is defined by bnear interpolation: U 


Pu 2 j — Uj, Pti 2j + 1 = ~(Wj + U J+l ) 


(2.3.3) 


Overbars indicate coarse grid quantities. A restriction 
following weighted average 


operator R : U — ► U is defined by the 


Ru i - 4 u 2j—l + ~U 2 j + ^«2j+i 


(2.3.4) 

where - i. defined to be sero outside G. Note that the matrices P and R are related by 
~ 2 ^ > but Inis property is not essential. ~ 

The fine grid equation (2.2.4) must be approximated by a coarse grid equation 

Au = f 
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the Mowing. The fine grid problem (2.2.4) is equivalent to 


(Au,v) = (f,v), ueU,Vv&U 


(2.3.5) 

( \ the standard inner product on U. We want to find an approximated solution Pu 

Si! 

operator that may be different from P. 

(2.3.6) 


(APu,Pv) = (f,Pv), ueu,*veu 


0T (P'APu,i) = (F“/,®), ueU,Vv€V ( 2 - 3 - 7 ) 

where now of course (., .) is over V, and soperscript * denotes the adjoint (or transpose ,n 
this case). Equation (2.3.7) is equivalent to 


Au = f 


(2.3.8) 


with A = RAP < 2 ' 3 ' 9 > 

and / = R/; we have replaced F' by R. This choice of A is called Gofer-tin coorse 9 rid 

WUhTetn’d R given by (2.2.3), (2.3.3) and (2.3.4), Equation (2.3.9) results in the follow, ng 


Aui = h' 2 ( 2 u 1 - fi 2 ) 

Au-j = h~ 2 {-Uj-i +2uj - uj+i) , j = 2,3, 1 

Afi n = fi~ 2 (— fin-i + Un) 


(2.3.10) 


which is the coarse grid equivalent of the left-hand side of (2.2.3). Hence, in the present 
advantages, as we shall see. 

Coarse grid correction = A _ u is to be approximated 

Let u be an approximation to the solution of (2.2.4). 1 he error e-u a 

on the coarse grid. We have 12 3 111 

Ae = -r = Au - J v ' ' ’ 
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The coarse grid approximation u of -e satisfies 

Aii = Rr (2.3.12) 

In a two-grid method it is assumed that (2.3.12) is solved exactly, the coarse grid correction 
to be added to it is Pu : 

u:=u + Pu (2.3.13) 

Linear two-grid algorithm 

The two-grid algorithm for Unear problems consists of smoothing on the fine grid, approxima- 
tion of the required correction on the coarse grid, prolongation of the coarse grid correction 

to the fine grid, and again smoothing on the fine grid. The precise definition of the two-grid 
algorithm is 

comment Two-grid algorithm; 

Initialize u ° ; 

for i := 1 step 1 until ntg do 
u'/ 3 := S(u°,A,f, i/j); 
r :~ f — At* 1 / 3 ; 
u := A 1 Rr; 

u 2 / 3 := 1^/3 + p 

U 1 := S(uV 3 ,A,f,v 2 ); 

u° := u 1 ; 


The number of two-grid iterations carried out is ntg. S(u°, A, /,i/ a ) stands for ^ smoothing 
iterations, for example with the Gauss-Seidel method discussed earfier, appUed to Au = /, 
starting with u°. The first appUcation of S is called pre- smoothing, the second post- smoothing. 

Exercise 2.3.1 Derive (2.3.10) (Hint. It is easy to write down RAui in the interior and at 
the boundaries. Next, one replaces u, by Pu,.) 


2.4 Two-grid analysis 

The purpose of two-grid analysis (as of multigrid analysis) is to show that the rate of conver- 
gence is independent of the mesh-size h. We wiU analyse algorithm (2.3.14) for the special 
case V\ = 0 (no-presmoothing). 


Coarse grid correction 

From (2.3.14) it foUows that after coarse grid correction the error e 2 / 3 = u 2 > 3 - u satisfies 


e 2 / 3 _ e l/3 + p^l/3 _ p e l/3 


(2.4.1) 
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with the iteration matrix or error amplification matrix E defined by 


E = I-PA X RA 

(2.4.2) 

We will express e« 3 explicitly in terms of «>/» This is possible only in the present simple 
one-dimensional case, which is our main motivation for studying this case. Let 

1 /3 

e 1 / 3 = d + Pe, with ej = e 2 ' 

(2.4.3) 

Then it follows that ^ = ^ = £<J 

(2.4.4) 

We find from (2.4.3) that 


, 1 i/3 . J/ 3 _ L 1 / 3 

d,2j = 0, d2j+i — — 2 Cj e 2i+ 1 2 C2,+2 

(2.4.5) 

Furthermore, from (2.4.5) it follows that 


RAd = 0 

(2.4.6) 

so that ^ 2/3 _ 

(2.4.7) 


N^t, we consider the effect of post-smoothing by one Gauss-Seidel iteration. From (2.2.5) it 
follows that the error after post-smoothing e 1 = u - u is related to e y 


2 c} = eT 


-e}_ 1 +2e} = e% , j = 2,3, ...,2n- 1 

“4i-l + e 2n = ® 

Using (2.4.5)(2.4.7) this can be rewritten as 

A = 0 

e\j = 2^ 2j+1 + 4 e ^ J 2 ’ e *j +1 = 2^ ’ 

^2n = ^2n— 1 


(2.4.8) 


(2.4.9) 


By induction it is easy to see that 

=,1 . 

' 2 f - 3 


41 < llldlloo , MU = majcfldil : 3 = 1 » 2 »-» 2n > 


(2.4.10) 
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Since d — e 2 / 3 , we see that Ganss-Seidel reduces the maximum norm of the error by a factor 
2/3 or less. 

Rate of convergence 

1 /3 

Since e 0 =0 because of the boundary conditions, it follows from (2.4.5) that 

MU < “lAe 0 !^ (2.4.11) 

since e 1 ^ 3 = e° (no pre-smoothing). 

From (2.4.9) it follows that 



Ae \j 

- 3 ri 

— 4 1 “ 

16^2j-l 64^2j— 3 — 

(2.4.12) 


Ae \j-\ 

- 

~ 2 € 2j 

Hence, using (2.4.10), 

\ Ae 2j\ 

< jMI« . 


(2.4.13) 

Substitution of (2.4.11) 

gives 

Ik 1 Ik. < 


(2.4.14) 


where r = Ae is the residual. This shows that the maximum norm is reduced by a factor of 
5/12 or better, independent of the mesh-size. 


This type of analysis is restricted to the particular case at hand. More general cases will be 
treated later by Fourier-analysis. There a drawback is of course the assumption of periodic 
boundary conditions. The general proofs of rate of convergence referred to in the introduction 
do not give sharp estimates. Therefore the sharper estimates obtainable by Fourier analysis 
are more useful for debugging codes. On the sharpness of rate of convergence predictions 
based on Fourier analysis, see [24]. 

Again: the essential principle 

How is the essential principle of multigrid, discussed in Section 2.2, recognized in the foregoing 
analysis? Equations (2.4.6) and (2.4.7) show that 

RAe 2/3 = 0 (2.4.15) 

Application of R means taking a local average with positive weights; thus (2.4.15) implies 
that Ae 2 / 3 has many sign changes, and is therefore rough. Since Ae 2 / 3 = Au 2 / 2 - / is 
the residual, we see that after coarse grid correction the residual is rough. The smoother 
is efficient in reducing this non-smooth residual further, which explains the h-independent 
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reduction shown in (2.4.14). 

Exercise 2.4.1 In the definition of G (2.2.2) and G (2.3.1) we have not included the point 
x = 0, where a Dirichlet condition holds. If Neumann condition is given at x = 0, the point 
x = o’ must be included in G and G. If one wants to write a general multigrid program for 
both cases, x = 0 has to be included. Repeat the foregoing analysis of the two-grid algorithm 
with x = 0 included in G and G. Note that including * = 0 makes A non-symmetric. This 
difficulty does not occur with cell-centered discretization, to be discussed in the next chapter. 


3 Basic Iterative Methods 

3.1 Introduction 

Smoothing methods in multigrid algorithms are usually taken from the class of basic iterative 
methods, to be defined below. This chapter presents an introduction to these methods. A 
more detailed account may be found in [141]. 


Basic iterative methods 

Suppose that discretization of the partial differential equation to be solved leads to the lol- 
lowing linear algebraic system 

Ay = b (3.1.1) 


Let the matrix A be split as 

A = M - N (3-1.2) 

with M non-singular. Then the following iteration method for the solution of (3.1.1) is called 


a basic iterative method : 

My m+i - Ny m + b 

(3.1.3) 

or 

y m+ 1 = Sy m + Tb 

(3.1.4) 

with 

S = M _1 N , T = M -1 

(3.1.5) 

so that we have 



y m+1 

= Sy m + M~ 1 b, S = M~ l N , N = M - A 

(3.1.6) 

The matrix S is called the iteration matrix of iteration method (3.1.6). 
Basic iterative method may be damped , by modifying (3.1.6) as follows 



y* = Sy m + M“'b 

y m+1 = uy* + (1 -u)y m 

(3.1.7) 
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By elimination of y* one obtains 


y m+ 1 = sy + wM-'fr 


with 


(3.1.8) 


S -u>S+(l-u>)I (3.1.9) 

rel e ted S by ValUeS “ f ‘ he U " daraPed itera,i ° n matrix S and ,he da ">P«< Nation matrix S' are 

A(S ) = wA(S') + 1 - u> (3.1.10) 

Although the possibility that a divergent method (3.1.6) or (3.1.8) is a good smoother (a con- 

method^arftoTe f J M ^ ******' ^ m ° St Ukely Candidates for S ood s ™°thing 

methods are to be found among convergent methods. In the next section, therefore some 
m L“ rgenCG ° f baSiC iteratiVe meth ° dS are Presented - F ° r “- background see 

Exercise 3.1.1 Show that (3.1.8) corresponds to the splitting 

M* = M/u, N* = A-M* (3.1.11) 

3.2 Convergence of basic iterative methods 

Convergence 

- C ZZTT f ° r <3 ' 1 ' ^ ,hC M ° Wi " S C °"“ P ‘ S Play a ” * m P ortant role. We have 

my - iSy + b, so that the error e m = y m - y satisfies 


» m+1 = Se m 


(3.2.1) 


(3.2.2) 


As shown in many textbooks, we have 
Theorem 3.2.2 Convergence of (3.1.3) is equivalent to 

p(S < 1 

with p(S ) the spectral radius of S. 

Regular splittings and M- and ItT-matrices 

Definition 3.2.2 The splitting (3.1 .2) is railed regular if M- > 0 and N > 0 (elementwise) 
The splitting is convergent when (3.1.3) converges. ^ ' 

Definition 3.2.3 ([129], Definition 3.3). The matrix at is called an M-matri* if a ( , < 0 for 
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all ij with * ^ j, A is non-singular and A 1 > 0 (elementwise). 

Theorem 3.2.3 A regular splitting of am M-matrix is convergent. 

□ 

Proof. See [129] Theorem 3.13. 

tion (see later). V j convergent iterative method can always be turned into a 

smoothing Ztholto “Z meth ° dS 

to be discussed are obtained easily, using equations (3.1.8), (3.1.9) and (3. • ) 

Hon rp it is worthwhile to try to discretize in such way that the resulting matrix A is 
„ I. easy to see if a discretisation matrix is an M-ma,n* we 

present the following theorem. 

Definition 3.2.4 A matrix A is called irreducible if from (3.1.1) one cannot extract a sub- 
system that can be solved independently. 


Definition 3.2.5 A matrix A is called a K-matnx if 

> 0, Vi , 


a; 


dt 


hj 


and 


< 0, Vi, j with i^j 
Y dij > 0 , Vi , 


(3.2.3) 

(3.2.4) 

(3.2.5) 


□ 


with strict inequality for at least one i. 

Theorem 3.2.4 An irreducible A'-matrix is an M- matrix. 

Proof. See [141]. 

Note that inspection of the K - matrix property is easy 

The following theorem is helpful in the construction of regular splittings. 

Theorem 3.2.5 Let A be an M- matrix. If M is obtained by replacing certain elements «* 
by values k, satisfying «, < k, < 0, then A = M - iV .s a regular spbttmg. 
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Proof. This theorem is 
2.2 in [87]. 


an easy generalization of Theorem 3.14 in [129] suggested by Theorem 

□ 


The basic iterative methods to be considered all result in regular splitting anH t 
numerically stable algorithms, if A is an M-matrix. This is one reason why it is advisable 

maWxTs an Mm * trix e<, ” ati °" ‘° bC S ° lved in S “ ch a way * ha * resulting 

matrix is an M-matnx. This may require upwind differencing for first derivatives Another 

reason ,s the excluston of numerical wiggles in the computed solution, because a monotonicit, 

principle is associated with the M-matrix property. y 


3.3 


Examples of basic iterative methods: Jacobi and Gauss-Seidel 


tag spbutags 0 ^' m0,tly, basic iterative "'e«hods by defining the correspond- 

s^rulred grif A fr ° m a on a two-dimensional 


Point Jacobi. M = diag (A). 

H 4 A ° btained fr ° m A by replacin & for aU i,j With j?i,i± 1 by zero 
the forward ordering of the grid points of Figure 3.3.1 this gives horizontal line Jacobi- 

horizon^ n™", 'V ^ ° f FigUre 332 ° ne obtains vertical Une Jacobi One 

Jacobi lt6ratl0n f ° UOWed ^ ° ne V6rtiCal Une JaC ° bi iteratin ^es alternating 
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Figure 3.3.1: 
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Grid point orderings for point Gauss-Seidel. 


Point Gauss-Seidel. M is obtained from 


A replacing a tJ for all i,j with j > i by zero 
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Figure 3.3.2: Grid point orderings for block Gauss-Seidel 
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Block Gauss-Seidel. M 


is obtained from A by replacing a, } for all t, j with j > i + 1 by 


zero. 


From Theorem 4.2.8 it is immediately dear that, if A is an M-matrix, then the Jacobi and 
Gauss-seidel methods correspond to regular splittigs. 

orderings. q fficeg therefore , to discuss orderings of computational grid points. 

computational grid. 1 1 .1 r w h} c h j s enough to illustrate the basic 

We restrict ourselves to a two-dimensional grid G, whicti is enougn io 

ideas. G is defined by 

G = {(>, j) : i = 1,2 1,2, ( 3 ' 31 * 

The points of G represent either vertices or cell centres, depending on the discretization 
method. 

Forward or lexicographic ordering 

The grid points are numbered as follows 


b — 1 4- f 7 — 1 ^ / 


(3.3.2) 


Backward ordering 

This ordering corresponds to the enumeration 

jfe = IJ + 1 - i - (j ~ IK 
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White-black ordering 

This ordering corresponds to a chessboard colouring of G, 
and then the white points, or vice versa; cf. Figure 3.3.1. 


numbering first the black points 


Diagonal ordering 

The points are numbered per diagonal, starting in a corner; see Figure 3.3.1. Diferent variants 
are obtained by starting in different corners. If the matrix A corresponds to a discrete oper- 
ator with a stencil as in Figure 3.3.3(b), then point Gauss-Seidel with the diagonal ordering 
ol figure 3.3.1 is mathematically equivalent to forward Gauss-Seidel. 


(a) 


(b) 


(c) 


Figure 3.3.3: Discretization stencils. 


Point Gauss-Seidel- Jacobi 

We propose this variant in order to facilitate vectorized and parallel computing; more on this 
hortly. M is obtained from A by replacing a tJ by zero except a„ and . We call this 
point Gauss- Seidel- Jacobi because this is a compromise between the point Gauss-Seidel and 
Jacobi methods discussed above. Four different methods are obtained with the following four 
orderings: the forward and backward orderings of Figure 3.3.1, the forward vertical line or- 
ermg of Figure 3.3 2, and this last ordering reversed. Applying these methods in succession 
results in four-direction point Gauss-Seidel-Jacobi. 

White-black line Gauss-Seidel 

This CM be seen as a mixture of lexicographic and white-black ordering. The concept is best 
llustrated with a few examples. With horizontal forward white-black Gauss-Seidel the grid 
points are visited horizontal line by horizontal fine in order of increasing j (forward), while 
per ne he grid points are numbered in white- black order, cf. Figure 3.3.1. The lines can also 
be taken m order of decreasing j, resulting in horizontal backward white-black Gaus-Seidel. 
Doing one after the other gives horizontal symmetric white-back Gauss-Seidel. Doing one 
after the other gives horizontal symmetric white-black Gauss-Seidel. The lines can also be 
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taken vertically; Figure 3.3.1 illustrates vertical backward white-black Gauss-Seidel. Com ^_ 
Lth„r"„d 6 ve,.icaa symmetric whte-black Gauss-Seidel gives a—, MU 
Gauss-Seidel. White-black Une Gauss-Seidel ordering has been propos [ 

corresponding to lines in the grid are updated simub 
uluS FoZard and Lctanf homontul line Gauss-Seidel correspond to the forward and 
£Ta!d' ordling, respectively, in Figure 3.3.1. Figure 3.3.2 gives some more order, ngs for 

block Gauss-Seidel. 

Symmetric horizontal line Gauss-Seidel is forward horizontal Une Gauss-Seidel foUowed 
by backward horizontal Une Gauss-Seidel, or vice yersa.d«erno<m 9 zebra ' “ 

horizontal zebra foUowed by vertical zebra Gauss-Seidel, or vice versa. Other combinations 

come to mind easily. 

S: wit “ttr^trdTscS above differ in their suitably for computing with 
vector or paraUel machines. Since the updated quantities are mutually independent, Jacob 

.hi points inside a diagonal are mutually independent if the structure^ of A 1 is £- . by 
Figure 3.3.3(b), if the diagonals are chosen as ,n Figure 3.3.1. The same “ ,ra * " ‘ J 

I Unis’ lithe grid* can be handled in parallel; for examp.e with 

^e fo^lrd ordering of Figure 3.3.1 the points on vertic* Hues Gauss-Sejde, pom.s o^ 
same colour can be updated simultaneously, resulting ,n a vector length of 1/2 or J/ 2, as 
case may be. 

Exercise 3.3.1 Let A = L + D + U, with [y = 0 for j > i, D = diag (A ), and «y = 0 for 
j > i. Show that the iteration matrix of symmetric point Gauss-Seidel is given by 

S = (U + D)-'L(L + D)- l U ( 3 - 3 ‘ 4) 


Exercise 3.3.2 Prove Theorem 3.3.1. 
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* 


Figure 3.3.4: Five-point stencil. 


3.4 


Exam ples of basic iterative methods: incomplete point LU factoriza- 


Complete LU factorization 

When solving Ay - b directly, a factorization A = LU is constructed, with L and U a lower 
and an upper triangular matrix. This we call complete factorization. When A represents a 
discrete operator with stencil structure, for example, a* in Figure 3.3.3, then L and U turn 

o e much less sparse than A, which renders this method inefficient for the class of 
problems under consideration. 

Incomplete point factorization 

W,tbmc°mpUte factorization or incomplete LU factorization (I LU) one generates a splitting 
L IpdU: " g SParSe a " d ‘° COmPUte IOWCT Md Upper ‘«-ianguIar factors 

M = LU 

If A is symmetric one chooses a symmetric factorization: 


An alternative factorization of M 


is 


M = LL t 


M = LD l U 


(3.4.1) 

(3.4.2) 


(3.4.3) 

7ianmT- P n t€ P °iT fa *° riz f° n ' D is chosen to ^ a diagonal matrix, and diag ( L ) = 

follows A arwhC 3 T d ( n‘ 4J) ^ equiva]ent ‘ L D and U ar « determined as 

follows A graph Q of the incomplete decomposition is defined, consisting of two-tuples (i j) 

^or w ic tee ements l t] , d u and u tJ ae allowed to be non-zero. Then L, D and U are defined 

{LD~'U) kl = a kl , V(Ar, /) e Q ( 3 . 4 . 4 ) 

We will diseuss a few variants of ILU factorization. These result in a splitting A - M — N 
‘ incomplete point factorization is obtained if D is defined by 
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, . . v. , f n J- rt h with a G 2R a parameter, and D a diagonal matrix defined 

by 4 il = S From now o'n the modified version will be discussed, since the unmodified 

• r if fc c a CTwia l rase This or similar modifications have been investigated in the 

ZZ oSugrid methods in‘[65], [97], 183], [82] aud [145], [147], A srvey is given in [142], 
We will discuss a few variants of modified ILU factorization. 

LllThrgHdlfgiven by (3.3.1), let the grid points be ordered according to (3.3.2), and let 
the structure of the stencil be given by Figure 3.3.3. Then the grap o 

Q = {(k,k — I), ( k,k - 1), (M). {k,k+ 1), (k,k+ /)} (3.4.5) 

For brevity the Mowing notation is introduced 

a k = a kk -i, c k = o fe , fc _i, d k = a kk , q k = a kM u 9k = a kM i 

Let the graph of the incomplete factorisation be 

f b V* “S,.' "Because tile graph contains live elements, the resulting 

mertod is called ’’file-point ILU. Let be the IJ , IJ matrices with elements m. 

respectively, and similarly for u g. Then one can wrrte 

LD~'U = a + 7 + b + M + V + a ^”V + aS 1t ? + I 6 M + 7^ V 

From (3.4.4) it follows 


(3.4.6) 


a = a, 7 = c, M = <7i V — 9 
and, introducing modification as described above, 

6 + a6~ l g + c6~' g = d + ad. 

The rest matrix N is given by 

N = aS' 1 q + cti' 1 g + ad 
The only non-zero entries of N are 

n kk -I + 1 = n fc.fc+/-i = c ^fc_i5fc-i 

n kk = a(\n ktk -]+\\ + I) 

Here and in the Mowing elements in which indices outside the range [1, IJ] occur are to be 
replaced by zero. From (3.4.9) the Mowing recursion is obtained. 

(3.4.12) 


(3.4.7) 

(3.4.8) 

(3.4.9) 

(3.4.10) 

(3.4.11) 


6 k = d k - a k 8 k l,g k -i - c k b k ^q k - 1 + n kk 
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This factorization has been studied in [41]. 

I'T “ foUows that ° overwrite d, so that the only additional storage required 

.s for N. When requ.red, the res, dual b - Ay-*' can be computed as follows without using 

6 - Ay-*' = N ( »”+> - y) (3.4.13) 

which follows easily from (3.1.3). Since N is usually more sparse than A, (3.4.13) is a cheap 

way to compute the residual. For all methods of type (3.1.3) one needs to tore o ly M and 
I\ , and A can be overwritten. J 

Seven-point ILU 

The terminology seven-point IL f/ indicates that the graph of the incomplete factorization has 
seven elements. The graph Q is chosen as follows: 

6 = HM±/), (M±/=Fl), (Mil), (k ,k)} (3.4.14) 

For the computation of L,D and U see (141). L,DandU can overwrite A. The only 
needed. ^ re< "" " ^ ° ’’ if 0,16 prefers ’ elem “ ts ° f N «n computed when 

Nine-point ILU 

The principles are the same as for five- and seven-point ILU. Now the graph Q has nine 
elements, chosen as follows 6 P 

C = S,U{(M±/±i)} (3.4.15) 

with given by (3.4.14). 

For the computation of L, D and U see [141]. 

Alternating ILU 

Alternating IL U consists of one ILU iteration of the type just discussed or similar, followed 
by a second ILU iteration based on a different ordering of the grid points As an examnlc let 
the grid be defined by (3.3.1), and let the grid points be numbered actrdfng to 


k = IJ+\-j-[i- l )j 


(3.4.16) 


This ordering is illustrated in Figure 3.4.1, and will be called here the second backward or- 
ermg, o istinguish it from the backward ordering defined by (3.3.3). The ordering (3 4 16) 
turn out to be preferable m applications to be discussed later. The computation of the 

auTs P :L“i" faCt 7 a *r fac *°; s 1 • D a " d 0 is dis ™ ss ' d ta [Ml]. If alternating 

required for’i ’fl Id T n * “ in thc piace of A ’ 50 additional storage is 

required for L, D and U. N can be stored, or is easily computed, as one prefers. 
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Figure 3.4.1: Illustration of second backward ordering. 


OtTerTL^varian. are obtained for other choices of Q. See [88] for some posjbffita. In 

=2 — 

that the elements of A outside Q are subtracted from N. 

The following algorithm compotes an ILU factorisation for general £1 by incomplete Gauss 
elimination. A is an n x n matrix. We choose drag (L) - diag (17). 

Algorithm 1. Incomplete Gauss elimination 

A 0 := A 

for r := 1 step 1 until n do 
begin a r rT := sqrt (a^ 1 ) r-1 r 

for j > r A (r, j) G G do a r T j := /a rr 
for i > r A (t,r) € £? do ay r := a T ir /a T rr 
for ( i,j ) G G A t > r A j > r A (t,r) G G A(r,j) £ G do 

ah := a"" 1 - «5 


*1J U 

od od od 


tr^rj 


end of algorithm 1. 

Stokes equations in the vorticity-stream function formulation). 

Final factorizations and numerical stability of the associated algorithms has been 

E Tin 1871 if A is an M-matrix; it is also shown that the associated splitting is regular, so 
proved m [87] if A is an iM information on efficient implementations 
that ILU converges according to Theorem 4.2-4. and 

of ILU on vector and parallel computers, see [69], [68], [116], LlWj.lilO], l j, l J 
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Exercise 4.4.1 Derive algorithms to compute symmetric ILU factorizations A = LD 1 L T - 
N and A = LL — N for A symmetric. See [87]. 

Excise 4.4.2 Let A = L + D + V, with D = diag (A), k, = 0, j > i and Uij = 0, j < 
bhow that (3.4.3) results in symmetric point Gauss-Seidel (cf. Exercise 3.3.1). This shows 
that symmetric point Gauss-Seidel is a special instance of incomplete point factorization. 

3.5 Examples of basic iterative methods: incomplete block LU factoriza- 
tion 

Complete line LU factorization 

The basic idea of incomplete block LU- factorization (IBLU) (also called incomplete line LU- 
factonzation (ILLU) in the literature) is presented by means of the following example. Let 
the stencil of the difference equations to be solved be given by Figure 3.3.3(c). The grid point 
ordering is given by (3.3.2). Then the matrix A of the system to be solved is as follows: 

B i U x \ 

L*2 Z?2 U 2 

(3.5.1) 

Uj _ i 

Lj Bj j 

with Lj,Bj and Uj I x I tridiagonal matrices. 

First, we show that there is a matrix D such that 

A = (L + D)D~ 1 (D + U) (3.5.2) 

where 



L 


D 



(3.5.3) 
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We call (3.5.2 a line LU factorization of A, because the blocks in L,D and U correspond to 
(in our case horizontal) lines in the computational grid. From (3.5.2) it follows that 

A = L + D + U + LD~'U (3-5.4) 


One finds that LD~'U is the following block-diagonal matrix 

/ 0 ^ 

L 2 Dj l U 1 


LD~ X U = 


V 


LjDj^Uj J 


From (3.5.4) and (3.5.5) the following recursion to compute D is obtained 
D\ = B\, Dj = Bj - LjD-^Uj , j = 2,3, ...,«/ 


(3.5.5) 


(3.5.6) 


Provided Dj 1 exists, this shows that one can find D such that (3.5.2) holds. 


Nine-point IBLU , A 

The matrices Dj are full; therefore incomplete variants of (3.5.2) have been proposed. A 

incomplete variant is obtained by replacing LjD^Uj in (3.5.6) by its tridiagona part (i.e. 
replacing all elements with indices i,m with m ± i,i± 1 by zero): 


D x = Bi , Dj = Bj- tridiag {LjDjl x u i) 

(3.5.7) 

The IBLU factorization of A is defined as 


a = (L + b)b~\b + u) — n 

(3.5.8) 

There are three non-zero elements per row in L, D and (/; thus we call this 
For an algorithm to compute D and D see [141]. 

nine-point IBL U. 

The IBLU iterative method 

With IBLU, the basic iterative method (3.1.3) becomes 


r = b-Ay m 

(3.5.9) 

(L + D)b~\b + U)y m+ ' =r 

(3.5.10) 

y m+\ ._ y m+1 + y m 

(3.5.11) 

Equation (3.5.10) is solved as follows 


Solve (L + D)y m+l - r 

(3.5.12) 
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r := Dy m+1 
Solve (D + L)y n+l = r 

With the block partioning used before, and with Vj and r 3 denoting 
corresponding to block j , Equation (3.5.12) is solved as follows: 


(3.5.13) 

(3.5.14) 
/-dimensional vectors 


D,y^= ru j = 2, 3, ..., / 


(3.5.15) 


Equation (3.5.14) is solved in a similar fashion. 


Other IBLU variants 

Other IBLU variants are obtained by taking other graphs for L,D and U. When A corre- 
sponds to the five-point stencil of Figure 3.3.3, L and U are diagonal matrices, resulting in the 
five-pomt mLU variants. When A corresponds to the seven-point stencils of Figure 3.3 3(a) 
(b), L and U are bidiagonal, resulting in seven-point IBLU. There are also other possibilities’ 
to approximate LjD^Uj by a sparse matrix. See [6], [33], [7], [99], [107] for other versions 
o IBLU; the first three publications also give existence proofs for Dj if A is an M-matrix* this 
condition is slightly weakened in [99], Vectorization and parallelization aspects are disc’used 


Exercise 3.5.1 Derive an algorithm to compute a symmetric IBLU factorization A = 
+ D)D (D + L T ) — N for A symmetric. See [33], 


3.6 Some methods for non-M-matrices 

When non-self-adjoint partial differential equations are discretized it my happen that the 
resulting matrix A is not an M-matrix. This depends on the type of discretization and 
the values of the coefficients. Examples of other applications leading to non-M-matrix dis- 
cretizations are the biharmonic equation and the Stokes and Navier-Stokes equations of fluid 
dynamics. 

Defect correction 

Defect correction can be used when one has a second-order accurate discretization with a 
matrix A that is not an M-matrix, and a first-order discretization with a matrix B which is 
an M-matrix, for example because B is obtained with upwind discretization, or because B 
contains artificial viscosity. Then one can obtain second-order results as follows. 

Algorithm 1. Defect correction 
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begin Solve By — b 

for i := 1 step 1 until n do 
By = b - Ay 4- By 

y ■= y 

od 

end of algorithm 1. 


It suffices in practice to take n — 1 or 
already y has second-order accuracy, 
can be used to solve for y. 


2. For simple problems it can be shown that for n - 1 
B is an M- matrix; thus the methods discussed before 


Distributive iteration 

Instead of solving Ay = b one may also solve 


ABy = b, y = By 


(3.6.1) 


This may be called post- conditioning, in analogy with preconditioning, where one solves 
BAy = Bb. B is chosen such that AB is an M- matrix or a small perturbation of an 

M-matrix, such that the splitting . 

AB = M - N (3.6.2) 

leads to a convergent iteration method. From (3.6.2) follows the following splitting for the 

original matrix A A = MB~'-NB~' (3.6.3) 

This leads to the following iteration method 

MB- X y m+A = Nfl-'t/" 1 + b (3-6-4) 


y m-\ = y m + BM~ l {b - Ay m ) (3.6.5) 

The iteration method is based on (3.6.3) rather that on (3.6.2), because if M is modified so 
that (3.6.2) does not hold, then, obviously, (3.6.5) still converges to the right solution if it 
converges. Such modifications of M occur in applications of post-conditioned iteration to the 
Stokes and Navier-Stokes equations. 

Iteration method (3.6.4) is called distributive iteration, because the correction M (6 - 
Av m ) is distributed over the elements of y by the matrix B. A general treatment of this 
approach is given in [144], [146], [148], [150], [149], where it is shown that a number of well 
known iterative methods for the Stokes and Navier-Stokes equations can be interpreted as 
distributive iteration methods. 
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Taking B A and choosing (3.6.2) to be the Ganss-Seidel or Jacobi splitting results in 
the Kaczmarz [78] or Cimm.no [32] methods, respectively. These methods converge for every 

U nroof of ’tb^ 7 and Jacobi “—ge for symmetric positive definite matrices 

a proof of this elementary result may be found in [70], Convergence is, however, usually 


4 Smoothing analysis 

4.1 Introduction 

The convergence behaviour of a multigrid algorithm depends strongly on he smoother. The 
efficiency of smoothing methods ,s problem-dependent. When a smoother is efficient for a 

\Z l rl7 V T 'l, CaUed r ° bUSL ™ S C ° nCe|>t wiU be "*<•« "X>re precise shortly 
a certain class of problems. Not every convergent method has the smoothing property 

but fo, symmetric matrices it can be shown that by the introduction of suitable amount of 

damping every convergent method acquires the smoothing property. This property says little 

about the actual efficiency A convenient tool for the study of smoothing efficiency is Fourier 

is th ’ * , aPPbe ' 1 ‘° ‘ he """-Vminetrir case. Fourier smoothing analysis 

IS the mam topic of this chapter. 6 y 

Many different smoothing methods are employed by users of multigrid methods. Of course 
m order to explain the basic principles of smoothing analysis it suffices to discuss only a few 

method'f ^ Way t ° f ! IUStra ;.' 0n - T ° fac,litate the making of a good choice of a smoothing 
method for a particular application it is, however, useful to gather smoothing analysis results 

wwt ^ n'T thr ° Ugh the htera.ture in one place, and to complete the information 
where results for important cases are lacking. 

4.2 The smoothing property 

The smoothing method is assumed to be a basic iterative method as defined by (3 1 3) We 

Theore nVv! A “ t- A ' n ’ atriX , 0f,en ' tbe smoother is obtaiMd « b " described in 
lheorem 3.2.5, in practice one rarely encounters anything else. 

The smoothing property is defined as follows ([57]): 

4 ' 2 ; 1 Smoothing property. S has the smoothing property if there exist a 
constant C s and a function V (u) independent of the mesh-size h such that 


HASH < C s h-* m r,(v), t](u) 0 for 

where 2m is the order of the partial differential equation to be solved. 


(4.2.1) 
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Here 5 is the iteration matrix defined by (3.1.5). The smoothing property implies converse 
[141] but, as already remarked, the converse is not true. In [145] it is shown tha a conve g 
method can be turned into a smoother by damping; for a fuller discussion see [1 ]. 

fn SrthTsmoothing property is shown for a number of iterative methods The smooth- 
ing property of incomplete factorization methods is studied in [145], [147 . Non-symmetnc 
problems can be handled by perturbation arguments, as indicated by [57]. When the non- 
symmetric part is dominant, however, as in singular perturbation problems this doe no 
lead to useful results. Fourier smoothing analysis (which, however, also has its limitations) 
can handle the non-symmetric case easily, and also provides an easy way to optimize value 
of damping parameters and to predict smoothing efficiency. The introduction of damping 
does not necessarily give a robust smoother. The differential equation 

eter such that when it tends to a certain limit, smoothing efficiency deteriorates. Examples 
and' further discussion of robustness will follow. We will concentrate on Fourier smoothing 

analysis. 

4.3 Elements of Fourier analysis in grid-function space 

As preparation we start with the one-dimensional case. 

The one-dimensional case 

Theorem 4.3.1. Discrete Fourier transform. Let / = {0, 1,2, ...,n-l}. Every u . I -*• JR 
can be written as 


m+p T 

u = £ c fc *(0fc), *j(9 k ) = exp(ijOk), 0 k = 2t rfc/n, J € I 


(4.3.1) 


/e= — r 


where p = 0, m = (» - l)/2 for , odd and p = 1, m = n/2 - 1 for n even, and 


c k = n 1 ^u^j{-0k) 
j=o 


(4.3.2) 


The functions $(6) are called Fourier modes or Fourier components. For a proof of this 
elementary theorem see [141]. 

The multi-dimensional case 
Define 


V>j(0) = exp(ij0) 


(4.3.3) 
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with j e /, o 6 0, with 

I ~ — iiiiiii ■■■ijd), j a = 0, 1,2 ,...,n a — 1, a = 1,2, (4.3.4) 

0 = {8 : 0 = (0],0 2 , = 2nk a /n a , 

k a = -m a ,-m, a + l,...,m Q + p a , a = 1,2, ...,d) 

Furthermore, 0 ’ ^ = ( "° " ^ ** ° d<3 *** Pa = ly m “ = n “/ 2 ~ 1 for 

d 

J 6 = 

O '— 1 


(4.3.5) 
even. 

(4.3.6) 


be^written i' 3 ' 2 ' DiSCFete F ° Urier transfori « m d dimensions. Every u : / - JR can 


with 


C 0V](8) 

0ee 


i€/ a=l 


(4.3.7) 


(4.3.8) 


For a proof see [141]. 

The Fourier series (4.3 7) is appropriate for rf-dimensional vertex- or cell-centered grids with 

nw°hl C , i y COn<liti ?" s ' For ,he “« ° f Courier sine or cosine series to accommodate 
Diriclilet or Neumann conditions, see [141]. 

4.4 The Fourier smoothing factor 

Definition of the local mode smoothing factor 

Let the problem to be solved on grid G be denoted by 

Au = f (4.4.1) 

and let the smoothing method to be used be given by (3.1.6): 

u := Su + M -1 /, S = M - N = A (4.4.2) 

According to (3.2.1) the relation between the error before and after „ smoothing iterations is 

e 1 = S^e 0 

We now make the following assumption. 


(4.4.3) 


29 



Assumption (i). The operator S has a complete set of eigenfunctions or loco! modes denoted 
by ip(9), 9 6 0, with 0 some discrete index set. 


Hence 


s^m = a* ww 

with A(0) the eigenvalue belonging to ij>{9). we can write 

e fl = £<«*(*)’ a = 0 ’ 1 

6>ee 


(4.4.4) 


(4.4.5) 


a " d ° btain ej = A*(0)e? 

The eigenvalue *(*) is also caUed the amplification factor of the local mode *(«). 
tot! assume that among the eigenfunctions *(») we somehow distmgursh between smooth 
eigenfunctions (0 G 0.) and rough eigenfunctions ( 9 € 0 r ): 


0 = 0 s U0 r , 0 s n 0 r = 0 


(4.4.6) 


We now make the following definition. 

Definition 4.4.1. Local mode smoothing factor. The local mode smoothing factor p of 
the smoothing method (4.4.2) is defined by 

p = sup{|A(9)| : « € e,} (J ' 17i 

Hence, after a smoothings the amplitude of the rough components of the error are multiplied 
by a factor p v or smaller. 

rmtr^Xrlrannlysis a useful too, for examining the ,^.1 — 
ing methods we must be able to easily determine p, and to choose 0 S such that error 
e = xl>(9), 960, is weU reduced by coarse grid correction. This can be done if Assump 

(ii) is satisfied. 

Assumption (ii). The eigenfunctions ip{9) of S are harmonic functions. 

This assumption means that the series preceding (4.4.5) is a Fourier series. When this is so 
pi also called the Fourier smoothing factor. In the next section we will give conditions such 
teumption (ii) holds, and show how , is easily determined; bu, first we d.scuss the 

choice of 0 r • 
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Aliasing 

Consider the vertex-centered grid G given by (4.4.8) with n a even, and the corresponding 
coarse grid G defined by doubling the mesh-size: 


G = {x 

€ IR d : x = jh, j 

O'l j •* jd) y h — (/ij , /l2) /l^), 

3a 0,1,2,..., 7i a , h a — 1 j 7i a , a =■ 1,2,. 

:,d} 

(4.4.8) 

H 

II 

€ IR d : x = jh, j 

0*1 > J*2) ** • j 3d)i h — (^1, ~hd )? 

jot 0,1,2,..., 71 a , h a — 1 f 71 a , Oc = 1,2,. 

■;d} 

(4.4.9) 


with n Q - n a j 2. Let d - 1 , and assume that the eigenfunctions of S on the fine grid G are 
the Fourier modes of Theorem 4.3.1: ^j(O) = exp (ij6), with 

0€0 = {0:0 = 27T*/n,, k = -m/2+1, -m/ 2 + 2, ...,m/2} (4.4.10) 

so that an arbitrary grid function v on G can be represented by the following Fourier series 

v j = J2 c ^( d ) (4.4.11) 

flee 

An arbitrary grid function v on G can be represented by 

(4.4.12) 

eee 

with i>(d) : G —> IR, ■ipj(O) = exp(ij6), and 

0 = {6:9 = 2nk/n 1 , k = -nj/2+ 1, -m/2 + 2, ..., m/2} (4.4.13) 

assuming for simplicity that m is even. The coarse grid point xj = jh coincides with the fine 
grid point x 2 j = 2jh. In these points the coarse grid Fourier mode i />(0) takes on the value 

= exp(ijO) = exp(i2j9) (4.4.14) 

For -n,/4+ 1 < k < m/4 the fine grid Fourier mode 4>(6 k ) takes on in the coarse grid points 
x^ the values of ^2 j{8k) = exp(2nijk/n l ) = V^Trfc/m), and we see that it coincides with 
the coarse grid mode ip(9 k ) in the coarse grid points. But this is also the case for another fine 
grid mode. Define k as follows 

0<k<n l /2: k' = - ni /2 + k 

-h 1 /2< k <0: k' = m/2 + k (4.4.15) 
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Then the fine grid Fourier mode tf(0 fc «) also coincides with 4>(0 k ) in the coarse grid points. 
On the coarse grid, 1>{0 k ') cannot be distinguished from VW- This is called aha ^^ 
rapidly varying function 0 k > ) takes on the appearance of the much smoother function VW 

on the coarse grid. 

Smooth and rough Fourier modes , 

Because on the coarse grid G the rapidly varying function cannot be approximated 

and cannot be distinguished from 4(9 k ), where is no hope that the part of the error w ic 
consists of Fourier modes *(*,), *' given by (4.4.15), can be approximated on the coarse 
grid G. This part of the error is called rough or non-smooth. The rough Fourier modes are 
defined to be t l>{0 k <), with k‘ given by (4.4.15), that is 

k ' € {-n a /2 + 1, + -n!/4}U{ni/4, m/4 + l,...,nj/2} (4.4.16) 

This gives us the set of rough wavenumbers 0 r = {0 : 9 = 2wk' /n, : k' according to (4.4.16)}, 
or 

0 r = {6 : 6 = 27rfc/ni, k = — «i/2 + 1, —n\j2 + 2, ..., ni/2 

and 6 € [— 7 T,— 7r/2]u[7r/2,7r]} (4.4.17) 

The set of smooth wavenumbers 0 S is defined as 0 a = 0\0 r? 0 given by (4.4.10) with d — 1, 
or 

0 S = {9 : 0 = 2tt fc/m, = -nj/2+1, -»i/2 + 2, »i/2 

and 0 G {—x/2, tt/2)} (4.4.18) 

The smooth and rough parts and of a grid function t> : G -> R can now be defined 
precisely by . £ ^(0), E .*(1) 

eee, ^ 9ee ’- (4.4.19) 

c* = E 

j - o 

In d dimensions the generalization of (4.4.17) and (4.4.18) (periodic boundary conditions) is 
Q = {6 :6 = (0 u 6 2 ,...,0d), 9 a = 2irk a /n a , k a = -n a /2 + 1, ...,n a /2} 


(4.4.20) 


0 3 = 0fl n (-ff/2, »r/2), 0r = © \ 0. 

a=l 
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Figure 4.4.1: Smooth (0 S ) and rough (0 r , hatched) wavenumber sets in two dimensions, 
standard coarsering. 


Figure 4.4.1 gives a graphical illustration of the smooth and rough wavenumber sets 
(4.4.20) for d = 2. 0 r and 0 S are discrete sets in the two concentric squares. As the 
mesh-size is decreased ( n a is increased) these discrete sets become more densely distributed. 

Semi-coarsening 

The above definition of 0 r and 0 S in two dimensions is appropriate for standard coarsening , 
i.e. G is obtained from G by doubling the mesh-size h a in all directions a = 1,2 
With semi-coarsering there is at least one direction in which h a in G is the same as in G. 
Of course, in this direction no aliasing occurs, and all Fourier modes on G in this direction 
can be resolved on < 5 , so hey are not included in 0 r . To give an example in two dimensions, 
assume h x = h x (semi-coarsering in the x 2 -direction). Then (4.4.20) is replaced by 

O s = © fi {[— 7r, 7r] x (— t/2, 7t/2)}, 0 a = 0 \ 0, (4.4.21) 

Figure 4.4.2 gives a graphical illustration. 

Mesh-size independent definition of smoothing factor 

We have a smoothing method on the grid G if uniformly in n a there exists a p * such that 

P<P'< 1, Vn a , a = 1,2, ...,d (4.4.22) 

However, p as defined by (4.4.7) depends on n a , because O r depends on n a . In order to 
obtain a mesh-independent condition which implies (4.4.23) we define a set 0 r D 0 r with 0 r 
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Figure 4.4.2: Smooth (0 4 ) and rough (0 r , hatched) wavenumber sets in two dimensions, 
semi- coarsening in a : 2 direction. 


independent of n a and define 


p = sup{|A(0)| : 6 G 0 r } 


(4.4.23) 


so that 


P < P 


(4.4.24) 


and we have a smoothing method if p < 1. For example, if 0 r is defined by (4.4.20), then we 
may define 0 r as follows: 


© r = n [-7T, tt] \ n c -»/2. */ 2 ) ( 4 - 4 - 25 ) 

a=l «= 1 


This type of Fourier analysis, and definition (4.4.23) of the smoothing factor, have been 
introduced by Brandt (1977). It may happen that A(0) still depends on the mesh-size, in 
which case p is not really independent of the mesh-size, of course. 

Modification of smoothing factor for Dirichlet boundary conditions 

If A (6) is smooth, then p - p = 0(/i™) for some m > 1. It may, however, happen that there 

is a parameter in the differential equation, say £, such that for example p- p - 0(h a /e). 

Then, for £ <C 1 (singular perturbation problems), for practical values of h a there may be 
a large difference between p and p. Even if p = 1, one may still have a good smoother. 

Large discrepancies between predictions based on p and practical observations may occut 

for singular perturbation problems when the boundary conditions are not periodic. It turns 
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out that discrepancies due to the fact that the boundary conditions are not of the assumed 
ype arise mainly from the presence or absence of wavenumber components 9 a = 0 (present 

ob serveTS f 8 7 HTl Tf'T' Dirichtet b ° U " dar * editions). It has been 

bserved [29], [83], [147] that when using the exponential Fourier series (4.3.7) for smoothing 

analysis of a practical case with Dirichlet boundary conditions, often better agreement with 

WaV6nUmberS ^ = ° Changing definiti ° n 

e D = { 0:9 = ( 9 l , 9 2 ,...,e d ), e a = 2 nk a /n a , k aJ t o, k a = -n a / 2 + i,...,» a /2} 

= 0 D n n(-/2, T/2), 0? = 0 D \ 0f (4.4.26) 

where the superscript D serves to indicate the case of Dirichlet boundary conditions The 
smoothing factor is now defined as 

PD = sup{|A(0)| : 9 <E 0? } (4.4.27) 

Figure 4.4.3 gives an illustration of Q r D , which is a discrete set within the hatched region for 
t7nex!7c7n SUPPOrt ^ ^ USefuln6SS ° f definitions (4-4.26) and (4.4.27) will be given in 
Notice that we have the following inequality 

PD <P< p (4.4.28) 

If we have a Neumann boundary condition at both = 0 and x a = 1, then 9 a = 0 cannot 
be excluded, but if one has for example Dirichlet at = 0 and Neumann at = 1 then the 
error cannot contain a constant mode in the direction, and 0 a = 0 can again be excluded. 

fa7u som' 4 ' 1 ! SU T 0Se " 1 = ^ Ch ' '' mesh - size of '■ mesh-size of G, one-dimensional 
case, p some integer), and assume periodic boundary conditions. Show that we have abasing 

9 k = 2 -xk/n u k 6 zn {(-n 1 /2,-n 1 /2fi]u[n 1 /2p,n 1 /2]} 

and define sets 0 r , 0 S . 

4.5 Fourier smoothing analysis 

Explicit expression for the amplification factor 

the Sm ° othin & factor P'P or PD according to definitions (4.4.7), (4.4.23) 

( ■ .27) we have to solve the eigenvalue problem St/>(9) = A (9)^(9) with S given by 
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Figure 4.4.3: Rough wavenumber set (0?, hatched) in two dimensions, with exclusion ot 
9 a = 0 modes; standard coarsening. 

(4,4.2). Hence, we have to solve WW = A(d)M««). In stencil notation (to be more fully 
discussed later) this becomes, in d dimensions, 


> eZ<l 


(4.5.1) 


with 7L — {0,±1,±2,...}. 

We now assume the following. 

Assumption (i). M(m, j) and N(m,j) do not depend on m. 

TM s assumption is satisfied ~ 
(ii) Of Section 4.4 is satisfied: the eigenfunctions of S are given y ( ■ • ). 

y N (j)exp{i(j + m)0] = exp(ime) £ N(j)exp(ij6) 

&2i d 

so that il>m{9) = exp(imO) satisfies (4.5.1) with 


A(0) = ^2 N{j)exp(ijO)/ E M(j)e x p(ij0) 


ieZ d 


ieZ d 


(4.5.2) 


p eri°didty requires that exp (im.S.) = exp[i(m„ + «„)*„], „ r exp{in j = Heil ce , £ 0 

as defined by (4.3.5) assuming to be even. Hence, the eigenfunctions are the Fourier 
modes of Theorem 4.3.2. 

Variable coefficients, robustness of smoother 

In general the coefficients of the partial differential equation to be solved will be variable of 
course. Hence Assumption (i) will not be satisfied. The assumption of uniform mesh-size is 
less demanding because ofen the computational grid G is a boundary fitted grid, obtained 
by a mapping from the physical space and is constructed such that G is rectangular and 
has uniform mesh size. This facilitates the implementation of the boundary conditions and 
of a multigrid code. For the purpose of Fourier smoothing analysis the coefficients M(m,j) 

^ *?• ^ ‘ frOZeT1 ’ We lna >" ex P ect to have a good smoother if p < 1 for all 

values M Oh N (j) that occur. This is supported by theoretical arguments advanced in [57l 
section 8.2.2. L J ’ 

A smoother is called robust if it works for a large class of problems. Robustness is a 

quantitative property, which can be defined more precisely once a set of suitable test problems 
has been denned. 

Test problems 

In order to investigate and compare efficiency and robustness of smoothing methods the 
following two special cases in two dimensions are useful 

- (ec 2 + s 2 )u iU - 2(e - l)csu 12 - (es 2 + c 2 )u 22 = 0 (4.5.3) 

~ £ ( u ,n + U 22 ) + cu 1 + su 2 = 0 (4.5.4) 

with c = cos/3, s = sin (3. There are two constant parameters to be varied: e > 0 and 3. 
Equation (4 5.3) is called the rotated anisotropic diffusion equation , because it is obtained by 
a rotation of the coordinate axes over an angle fi from the anisotropic diffusion equation: 


EU' U - u t 22 = $ 


(4.5.5) 


Equation (4 5.3) models not only anisotropic diffusion, but also variation of mesh aspect ratio 
ecause wit fi — 0,£ - 1 and mesh aspect ration fii//i 2 = 8~ l ! 2 discretization results in the 
same stencil as with £ = 6, h l /h 2 = 1 apart from a scale factor. With 3 ^ kw/2, k = 0 1 2 3 
(4 5.3) also brings in a mixed derivative, which may arise in practice because of the use of no’n- 
orthogonal boundary-fitted coordinates. Equation (4.5.4) is the convection-diffusion equation. 
It is not self-adjoint. For e « 1 it is a singular perturbation problem, and is almost hyper- 

bobc. Hyperbolic, almost hyperbolic and convection-dominated problems are common in fluid 
dynamics. 
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Equations (4.5.3) and (4.5.4) are not only useful for testing smoothing methods, but also 
for testing complete multigrid algorithms. General (as opposed to Fourier analysis) mult. gnd 
convergence theory is not uniform in the coefficients of the differential equation, and the the- 
oretical rate of convergence is not bounded away from 1 as £ | 0 or £ -> oo In the absence 
of theoretical justification, one has to resort to numerical experiments to validate a method, 
and equations (4.5.3) and (4.5.4) constitute a set of discriminating test problems. 

Finite difference discretization results in the following stencil for (4.5.3), assuming h x = 
h 2 = h and multiplying by h 2 : 


[A] = (£c 2 -M 2 )[- 1 2 -1] 
1 -1 

+ (e - 1 )cs 


— 1 
0 


2 

-1 


0 

-1 

1 


+ (ss 2 + c 2 ) 


-1 

2 

-1 


(4.5.6) 


The matrix corresponding to this stencil is not a ^-matrix (see Definition 3.2.6) if e- l)cs > 0. 
If that is the case one can replace the stencil for the mixed derivative by 


0 1 -1 
1 -2 1 

-1 1 0 


(4.5.7) 


We will not, however use (4.5.7) in what follows. 

A more symmetric stencil for [A] is obtained if the mixed derivative is approximated by 
the average of the stencil employed in (4.5.6) and (4.5.7), namely 


1 0 

0 0 

-1 0 


Note that for [A] in (4.5.6) to correspond to a K 
sc 2 + s 2 + (£ - l)cs > 0 and 


-1 

0 

1 

-matrix it is also necessary that 
es 2 + c 2 + (£ - l)cs > 0 


(4.5.8) 


(4.5.9) 


This condition will be 
With (4.5.8) there are 
have a AT-matrix. On 
the other two options 


violated if £ differs enough from 1 for certain values of c = cos f3, s = sin 
always (if (e- l)cs ^ 0) positive off-diagonal elements, so that we never 
the other hand, the ‘wrong’ elements are a factor 1/2 smaller than with 
Smoothing analysis will show which of these variants lend themselves 


most for multigrid solution methods. 


0. 
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Finite difference discretization results 
and multiplying by h 2 : 


in the following stencil for (4.5.4), with hj 


— fi2 = h. 


[A} = e 

-1 

-1 4 -1 

-1 

h r h 

+ C 2l-‘ 0 >) + *| 

r 1 

1 

O K-i 
1- 1 





(4.5.10) 


In (4.5.10) central differences have been used 
With upwind differences we obtain 


to discretize the convection terms in 


(4.5.4). 




-l 


4 

-1 


-1 


+ 2 [ ~ C ~ |c| C ~M] 


h 

+ 2 


M 

-s - js| 


(4.5.11) 


^fulfilled 10) glV6S a A matriX ° nly if the WeU known con ditions on the mesh Peclet numbers 

\c\h/e < 2, | s\h/e < 2 (4.5.12) 

wind'd-ff 4 ' 5 ' 1 ^ a n7 S - eSUltS ^ a /i '- matrix > whic h ^ the main motivation for using up- 
ZiTZ C ? ° f n ’. ln applications (for example, fluid dynamics) conditions (4 5 12) are 

c feat: M ^ iS hard ‘° hMdle With ■"'•hods; therefore dis- 

cretization (4.5.11) will mainly be considered. 


Definition of robustness 

We can now define robustness more precisely: 
above test problems, p < p* < 1 or < p* 

ho > h > 0. 


a smoothing methods is called robust if, for the 
< 1 with p * independent of e and h, for some 
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Numerical calculation of ^'^“‘^f^fflcal. to compute |A(.)|, and to find 
Using the explicit express! ( . • • ) ^ beMe , he Fourier smoothing factors p or 

its largest value on the discr r r , . 2 q) G r (4 4.21) various values of 

p D . By choosing in the definition of of the mesh- 

no one may gather numenca evi ence difficult numerically, since this 

independent smoothing factor p defined m 4 ^ can be found analytically, as 

involves finding a maximum on iuj^ £'• J 8|A(»)|/W, = 0,a = 1,2, ...,d 

- thaU see shortly. Extrem ^ ippUcat i 0 „ one can compute , for the 

and at the boundary of O r . Ut course, ror a p yy Bmit n _> qo. In the 

values of n a occurring in this app cation, wi ou^ _ ^ ig found that the smoothing 
Mowing, we often pmsen «s»Us /or a, ^ ^ „ those cases where , 

apprtc'S 8 An analysis wifi he given of what happens in those cases. 

All smoothing methods to be discussed in this chapter have been defined in Section 3.3 
to 3.5. 

Local smoothing _ wWp the coefficients are not 

Local freezing of the coefficients is not obtained as a boundary 

smooth. Such points ®y 0 « '^e compu^ ^ non . smoo th boundary. Near points on 
fitted coordinate mapping p y where the physical domain boundary is not 

the boundary which are the images o p nothine Derformance often deteriorates. 

points in a neighbourhood of t e “ s '" 6 “ " ^derations of vec ,or and parallel computing. 

W «■ analysed theoretically in ,110, and [34,. 

4.6 Jacobi smoothing 
Anisotropic diffusion equation 

Point^ Jacobi with damping corresponds to the fofiowing spotting (cf. Exercise 3.1..), in 
stencil notation: M(0) = „-M(0), M (y) = 0, jt 0 («-D 

Assuming periodic boundary conditions we obtain, using (4.5.9) and (4.5.2), in the special 
case c = 1, s = 0 = , + u(£ cos ,, _ £ + co s f» 2 - 1)/(1 + 0 ( 4 - 6 ' 2 > 
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Because of symmetry 0 r can be confined to the hatched region of Figure 4.6.1. Clearly, 
P > l/ > ( 7r ’ ^ ) I = U — 2a/| > 1 for u ^ (0, 1). For u> E (0, 1) we have for 8 E CDEF : \(w, n) < 



Figure 4.6.1: Rough wavenumbers for damped Jacobi. 

A(0, tt/2), or 1 — 2 ui < X(8) < 1 — a?/( 1 + £ ). For 6 E ABCG we have 

A(7t,7t/2) < A (8) < A(rr/2, 0), 

or 1 - [(1 + 2e)/(l + e))u < A ( 8 ) < 1 - [ £ /( 1 + £ )] w . 

Hence 

p = max{|l - 2w|, |1 - ^|, |1 - |1 - y-^l} (4.6.3) 

p = (2 + £ )/(2 + 3 £ ), w = (2 + 2 £ )/(2 + 3e) (4.6.4) 

For £ = 1 (Laplace’s equation) we have p = 3/5, w = 4/5. For £ < 1 this is not a good 
smoother, since bin p — The case e > 1 follows from the case e < 1 by replasing e by l/ £ . 

Note that p is attained for 8 E 0 r , so that here 

P = P (4.6.5) 

For o» = 1 we have p = 1, so that we have an example of a convergent method which is not a 
smoother. 

Dirichlet boundary conditions 

In the case of point Jacobi smoothing the Fourier sine series is applicable (see [141]), so 
that Dirichlet boundary conditions can be handled exactly. It is found that with the sine 
series A(0) is still given by (4.6.2), so all that needs to be done is to replace 0 r by 0(? in the 
preceding analysis. This is an example where our heuristic definition of p D leads to the correct 
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result. Assume », = n 2 = ». The whole of 0? is within the hatched region of Figure 4.6.1. 
Reasoning as before we obtain, for 0 < £ < 1: 

\(n,n)< A(0)< A( 27 r/n, tt/2), A(tt, tt/2) < A(0) < A(tt/ 2, 27r/n) (4.6.6) 

Hence p p = max{|l - 2u,|, |1 - £0,(1 + 1* 2 /n>)/(\ + e)\, so that p D = p + 0(n~ *), and again 
we conclude that point Jacobi is not a robust smoother for the anisotropic diffusion equation. 


Line Jacobi . _ T , . 

We start again with some analytical considerations. Damped vertical lme Jacobi iteration 

applied to the discretized anisotropic diffusion equation (4.5.6) with c = 1, s - 0 correspon s 
to the splitting r 

- 1 1 0 2 + 2e 0 | (4.6.7) 

-1 


[M] = u 


The amplification factor is given by 

A(0) = u>e cos 0i/(l + £ — cos O 2 ) + 1 — w (4.6.8) 

both for the exponential and the sine Fourier series. We note immediately that |A(tt, 0)| = 1 
if w = i, so that for u = 1 this seems to be a bad smoother. This is surprising, because as 
e | 0 the method becomes an exact solver. This apparent contradiction is resolved by taking 
boundary conditions into account. In Example 4.6.1 it is shown that 

p D = |A( 7 T,y>)| = £/(l + £ - cos <p) for u = 1 (4.6.9) 


where ip — 2?r / u. As 71 — ► oo we have 

pc ~ (1 + 27r 2 /i 2 /e) _1 (4-6.10) 

so that indeed lim p D = 0. Better smoothing performance may be obtained by varying u. In 
Example 4.6.1 it is shown that p is minimized by 

(4.6.11) 

3 + 2£ 

Note that for 0 < £ < 1 we have 2/3 < w < 4/5, so that the optimum value of u is only 
weakly dependent on e. We also find that for u in this range the smoothing factor depends 
only weakly on u. We will see shortly that fortunately this seems to be true for more general 

problems also. 

With u according to (4.6.11) we have 

p = (l+2£)/(l + 3£) (4.6.12) 
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Choosing w = 0.7 we obtain 


p = max{l - 0 . 7 /(l + £ ),0.6} (4.6.13) 

which shows that we have a good smoother for all 0 < e < 1, with an e-independent «. 

Example 4.6 1. Derivation of (4.6.9) and (4.6.11). Note that A (0) is real, and that we need 
consider ordy ^ ^ 0 It is found that dX/60, = 0 only for 0 X = 0,tt. Starting with p D , 
we see that max{|A(0)| : 0 € 0? } is attained on the boundary of 0^. Assume », = n 2 = n 
and define = 2* /». It is easily see that max{|A(0)| : 0 € 0?} will be either |A(<*,r/2)| or’ 
|A(7r,^)| If w - 1 it is |A(t, v»)|, which gives us (4.6.9). We will determine the optimum value 
tu no or p D ut for p. It is sufficient to look for the maximum of |A(0)| on the boundary 
ot 0 r . It is easily seen that J 

p = max{|A(0, 7 t/2)|, |A(7t, 0)|} = max{l - w/(l +£■), |l - 2u>\} 

which shows that we must take 0 < u < 1. We find that the optimal u> is given b, (4.6.11). 
Note that in this case we have p = p. y 

Equation (4.5.5), for which the proceeding analysis was done, corresponds to [3 = 0 in (4.5.3). 
or 0 x/2 damped vertical line Jacobi does not work, but damped horizontal line Jacobi 

S ou e use e general case may be handled by alternating Jacobi : vertical line followed 
y horizontal line Jacobi. Each step is damped separately with a fixed problem-independent 

cT 6 eXperimentation w = 0-7 was found to be suitable; (cf. (4.6.12) and 
(4.6.13). Table 4.6.1 presents results. Here and in the remainder of this chapter we take 

n, - n 2 - n, and 0 is sampled with intervals of 15°, unless stated otherwise. The worst case 
found is included in the tables that follow. 

Increasing n, or finer sampfing of 0 around 45° or 0°, does not result in larger values 
of p and than those listed in Table 4.6.1. It may be concluded that damped alternating 
Jacob, with a fixed damping parameter of u, = 0.7 is an efficient and robust smoother for the 
rotated anisotropic diffusion equation, provided the mixed derivative is discretized according 
to (4.5.8). Note the good vectorization and parallelization potential of this method 
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£ 

(4.5.6) 

(4.5.8) 

Pi PD 

a 

Pi PD 

P 

1 

0.28 

any 

0.28 

any 

10" 1 

0.63 

45° 

0.38 

45° 

10" 2 

0.95 

45° 

0.44 

45° 

10" 3 

1.00 

45° 

0.45 

45° 

10" 5 

1.00 

45° 

0.45 

45° 

10" 8 

1.00 

45° 

0.45 

45° 


Table 4 6 1- Fourier smoothing factors p,p D for the rotated anisotropic diffusion ^^uation 
“ d) discretized according to (4.5.6) or (4.5.8); damped aiternatmg Jacob, smooth, ng, 

u> = 0.7; n = 64. 

Convection-diffusion equation 

Forte'rveLion-diffnsion egua.ion discretized with stencil (4.5.11) the amplification factor 
of damped point Jacobi is given by 

m = c(2 cos 0, + 2 cos dr + ft."' + P^W + ft + ft) + 1 - " < 4 - 6 ' 14 > 

where Pi = ch/e. Pi = sh/c. Consider the special case: Pi = 0, Pr = 4/6. Then 

A(*,0)=l-u + W(l+‘) (4 ' 6 ' 15) 

so that |A(n,0)| - 1 as 6 1 0, for all to, hence there is no value of u, for which this smoother 
is robust for the convection-diffusion equation. 

Let us apply the line Jacobi variant which was found to be robust for t ^ e r ^ ted 
diffusion equation, namely damped alternating Jacob, w.,h a, - 0.7, 
diffusion test problem. Results are presented in Table 4.6.2. 

v f a rnmul ft — (1° and increasing n does not result in significant changes. 
Finer sampling of f) around fi - 0 ndmcw g be ^ dMnpcd 

Numerical experiments show u> - 0.7 to be a go _ , . b t an d 

easily, so that all in all is an attractive smoother. 
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e 

P 

0 

PD 

0 

1 

0.28 

0° 

0.28 

~w 

10- 1 

0.28 

0° 

0.29 

0° 

10- 2 

0.29 

0° 

0.29 

0° 

10- 3 

0.29 

0° 

0.29 

0° 

1(T 5 

0.40 

0° 

0.30 

0° 


Table 4.6.2: Fourier smoothing factors p,p D for the convection-diffusion equation discretized 
according to ( 4 . 5 . 11 ); damped alternating line Jacobi smoothing; w = 0.7; n = 64. 

Exercise 4.6.1 Assume semi-coarsening as discussed in Section 4.4: h, = h u h 2 = h 2 / 2. 
Show that damped point Jacobi is a good smoother for equation (4.5.5) with 0 < e < 1. 

Exercise 4.6.2 Show that bmp = 1 for alternating Jacobi with damping parameter w = 1 
applied to the convection-diffusion test problem. 


4,7 Gauss-Seidel smoothing 

Anisotropic diffusion equation 
Point Gauss-Seidel 

Forward point Gauss-Seidel iteration applied to (4.5.3) with c = 1, a = 0 corresponds to the 
splitting 


(4.7.1) 


0 


1 1 

-£ 2e + 2 0 

, [N] = 

Uj 

O 

O 

-1 


0 J 


[M] = 

The amplification factor is given by 

\(0) = ( £e ’ 92 + e te2 )/(-ee~ i$1 + 2e + 2 - e' 02 ) 
For s — 1 (Laplace’s equation) one obtains 


(4.7.2) 


P = |A(tt/ 2, cos 1 (4/5))| = 1/2 (4.7.3) 

To illustrate the technicalities that may be involved in determining p analytically, we give the 
details of the derivation of (4.7.3) in the following example. 


Example 4.7.1. Smoothing factor of forward point Gauss-Seidel for Laplace equa- 
tion. We can write 


l^(^)| 2 — (1 + cos /3)/(9 - 8 cos ^ cos — + cos j3) 

& 2 


(4.7.4) 
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with a = #1 + 0 2 , /3 = *1 - *2- Because of symmetry only a, P > 0 has to be considered. We 

have 9|A(0)|V0O = 0 for sin(a/2) cos(/3/2) = 0 (4 7 5) 

This gives a = 0 or a = 2* or /3 = *. For /J = rr we have a minimum: |A| ! = 0. With 
c = 0 we have |A W P = cos 2 W 2)/(2 - cos W 2)) 2 , which reaches a maximum for : /» = 2* 

1 e. at the boundary of 9,. With a = 2l we are also on the boundary of ft.. Hence, the 
maximum of |A(*)| is reached on the boundary of 0 r . We have W*/2,M ~ 

s i„«, _ 4 cos 0 2 ), of which the « 2 derivative equals 0 of 8 A cos « 2 4sm9 2 , 

0 2 = — ff/2, which gives a minimum, or 0 2 = ±cos (4/5). The largest maximum ,s ob- 
tained for 0 2 = cos-' (4/5). The extrema of |A(*,9 2 | are studied in similar fashion. Since 
aT«T« 2 ) = A(9 2 .9, ) there is not need to study |A(6„x/2)| and |A(9„x)|. Equation (4.7.3) 

follows. 

We will not determine p analyticaUy for £ # 1, because this is very cumbersome. To do this 
numerically is easy, of course. Note that Urn A(rr,0)= 1, Um A(x,0) = -1, so that forward 
point Gauss-Seidel is not a robust smoother for the anisotropic diffusion equation, if standar 
coarsening is used. See also Exercise 4.7.1. 

With semi-coarsening in the x 2 direction we obtain in Example 4.7.2: p < {( 1 + 0/(5 + 

e)Yl\ which is satisfactory for e < 1. For £ > 1 one should use semi-coarsenmg m the 
x, -direction. Since in practice one may ave £ < 1 in one part of the domain an £ 
another, semi-coarsening gives a robust method with this smoother only ‘V e ,Un “^dard 
semi-coarsening is varied in the domain, which results in more complicated code than standard 

multigrid. 

Example 4.7.2. Influence of semi-coarsening. We will show 

p < [(1 + 0/(5 + 0] 1/2 


(4.7.6) 


for the smoother defined by 94.7.1) with semi-coarsenmg m the r, Fron ' ( . 

i, follows that one may write |A(0)|- 2 = 1 + (2 + *)/<(*) with ^ =* < 2 + „ ha ' ve 

2 cos 0 2 )/[l + £ 2 + 2£ cos (0! - 6 2 ]. In this case, 0 r is given m Figure 4.4.2. On O r we have 

p(9) > (2 + 2£ - 2£cos Oi - 2 cos 0 2 )/(l + e) 2 > 2 /0 + £ ) ■ 

Hence |A(<9)| > [1 + 4/(1 + e)] -1/ 2 , and (4.7.6) follows. 

For backward Gauss-Seidel the amplification factor is A(-0), with X(9) given 80 

that the amplification factor of symmetric Gauss-Seidel is given by A(-0)A(0) From (4 7 2) 
it follows that |A(0)| = |A(-0)|, so that the smoothing factor is the square of the smoo g 
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factor for forward point Gauss-Seidel, hence, symmetric Gauss-Seidel is also not robust for 

the anisotropic diffusion equation. Also, point Gauss-Seidel-Jacobi (Section 3.3) does not 
work for this test problem. 


The general rule is: points that are strongly coupled must be updated simultaneously. Here 
we mean by strongly coupled points: points with large coefficients (absolute) in [Al. For 
example, in the case of Equation (4.5.5) with e < 1 points on the same vertical line are 
s rongly coupled. Updating these points simultaneously leads to the use of line Gauss-Seidel. 


Line Gauss-Seidel 

Forward vertical line Gauss-Seidel iteration applied to the anisotropic diffusion 
(4.5.5) corresponds to the splitting 


equation 


[M] 


-1 


o 1 

-£ 2e + 2 0 

, [N] = 

o 

o 

-1 


» J 


The amplification factor is given by 


(4.7.7) 


A(0) — ee t9i / ( 2e + 2-2 cos 0 2 — ee l9x ) (4.7.8) 

and we find Example 4.7.3, which follows shortly: 

p = max{5 -1 / 2 , ( 2/e + l) -1 } (4.7.9) 

Hence, Ujnp = 5" 1 / 2 . This is surprising, because for £ = 0 we have, with Dirichlet bound- 
ary conditions, uncoupled non-singular tridiagonal systems along vertical fines, so that the 
smoot er is an exact solver, just as in the case of line Jacobi smoothing, discussed before. The 
behaviour of this smoother in practice is better predicted by taking the influence of Dirichlet 
boundary conditions into account. We find in Example 4.7.3 below: 


£ < ' ( 1 T \/5)/2 : pd _ £[e 2 -f- (2s: + 2 — 2 cos </>) 2 ] J / 2 
€ > (1 + \/5)/2 : p D = e[e 2 + (2e + 2){2e + 2 — 2e cos y?)]' 1 / 2 (4.7.10) 


with <p - 2nh, h = 1/n, assuming for simplicity nj = n 2 = n. 
this can be approximated by 


For e < (1 + \/5)/2 and h J. 0 


PD ^ [1 + (2 + <p 2 /e) 2 } !/ 2 (4.7.11) 

and we see that the behaviour of p D as £ J 0, h | 0 depends on y> 2 / £ = 4 tt 2 h 2 /e For h I 0 
with £ fixed we have p D £ p and recover (4.7.9); for £ ] 0 with h fixed we obtain p D a 0. To 
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give a practical example, with h = 1/128 and £ = 1CT 6 we have p D 5! 0.0004. 

Example 4.7.3. Derivation of (4.7.9) and (4.7.10). It is convenient to work with 
|A(0)| -2 . We have 


|A(0)| -2 = [(2e + 2 - £ cos - 2 cos 0 2 ) 2 + e ‘ 


sin 


0i)/e 2 


Min JIA(0)|- 2 ' 6 € 0?} is determined as follows. We need to consider only 9 a > 0 It is 
found that dim-VMi = 0 for h = 0 for h = 0,x only. Hence the mm, mum ^U- 
on the boundary of 0? . Choose for simplicity n, = m = «, and define v - 2,/n. It easily 

seen that in 0^? we have 


iA(0i,^r 2 

| A ( tt , 0 2 )|- 2 
iA(7r/2,0 2 )r 2 


> 

> 

> 


| A ( tt / 2 , v )| 2 
|A(tt,¥5)| , 

| A ( 7 t / 2,^)|- 2 


, |A(vp,<? 2 )|- 2 > |A(^,7r/2)r 2 
|A(0!,7r/2)r 2 > |A(<^,7 t/2)| 2 , 
, |A(0i, 7r)| -2 > | A(</>, 7r)| — 2 


For £ < (1 + V5)/ the minimum is Wx/2,».)|- J ; for £ > 1 + t/5)A 

nL jr/2)| -3 . This gives ns (4.7.10). We continue with (4.7.9). The behaviour of |A(«)| on 
the boundary of 0. is found simply by letting », - 0 I in the preceding results. No. there , 
also the possibility of a minimum in the interior of 0. because H = 0 .. flowed, but 
leads to the minimum in (x/2,0), which is on the boundary, and (4.7.9) foUows. 

Foliations (4 7 9) and (4.7.10) predict bad smoothing when e > 1. Of course, for £ » 1 
horizontal Lie Gauss-Seidel should be used. A good smoother for arbitrary £ is alternahng 
line Gauss-Seidel. For analytical results, see [141). Table 4.7 1 presents numerical values of 
p and p D for a number of cases. We take ni = n 2 = n = 64, 0 - 1st / 12, k - 0, 1, . •••» 
(4.5.3), and present results only for a value of /? for which the largest p or pr> is o aine . 
the cases listed, p = PD- Alternating line Gauss-Seidel is found to be a robust smoother for 
the rotated anisotropic diffusion equation if the mixed derivative is discretized according o 
(4.5.8), but not if (4.5.6) is used. Using under-relaxation does not change this conclusion. 

Convection-diffusion equation 

Point Gauss-Seidel . r 

Forward point Gauss-Seidel iteration applied to the central discretization of the convec 

diffusion equation (4.5.10) is not a good smoother, see [141]. 

For the upwind discretization (4.5.11) one obtains, assuming c > 0, s > 0: 


m = 


e'Ml +(| fi| - Pi)/2] + e ,fl2 [l + (lF 2 l-P 2 )/2] (4.7.12) 

4 + |Pi| + |P 2 | - e-^ 1 [1 + (Pi + l^il)/ 2 ! - e'^l 1 + ( P 2 + |f > 2|)/2] 

with Pi - ch Is P 2 = sh/e the mesh-Peclet numbers (for simpUcity we assume m =n 2 ). 

Fm pfVo ?\< o we hive W 0, x)| = |ft/(« - ftM, •»* «o 1 as |P 2 | - oo. To avo.d 
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£ 

(4.5.6) 

(4.5.8) 

Pi Pd 

P 

PiPD 

P 

1 

0.15 

any 

0.15 

any 

10- 1 

0.38 

105° 

0.37 

15° 

10- 2 

0.86 

45° 

0.54 

15° 

10~ 3 

0.98 

45° 

0.58 

15° 

1(T S 

1.00 

45° 

0.59 

15° 


V’ 1: f F ° U , ner sm ° othln g factors p,p D for the rotated anisotropic diffusion equation 
M.3) discretized according to (4.5.6) and (4.5.8); alternating line Gauss-Seidel smoothing; 


his the order m which the grid points are visited has to be reversed: backward Gauss-Seidel 
Symmetric point Gauss-Seidel (forward followed by backward) therefore is more promising 
for ^he convection-diffusion equation. Table 4.7.2 gives some numerical results for o for 

n \~ n \ = , 64 ' We g ' Ve FeSultS f ° r a Value of P in the {P = Aj7t / 1 2 : k = 0 1 2 ... 23) for 
which the largest p and pjj are obtained. ? ’ ’ 

Although this is not obvious from Table 4.7.2, the type of boundary condition may make 

7" e : ence ; Fo ; n ins ; aace ’ for £ = 0 and * i ° ™e finds numerically for forward point 
SS ei e • P | (0,7r/2)| - 1/ \/5, whereas Ump D = 0, which is more realistic, since as 

bv ° I 116 rrr hCT u beCOmeS “ 6XaCt SOlVer ' Th ' e differenCe betWeen P and PD is explained 
7= 2^6 1 = V = 2 * h and < 1 we have |A(*>, ,r/2)| 2 = 1/(5 + y + ly 2 ) with 

point Calis^S 7M 105 ° X abIe 4 7 2 Sh ° WS rather large smoothing factors. In fact, symmetric 
p int Gauss-Seidel smoothing is not robust for this test problem. This can be seen as follows, 
li r\ < U, > 0 we find 


A(f,0) 


1 ~ P\ + i 1 + P 2 - i 

3 — Pi + f 3 - Pj + P2 - * ( 1 - Pi ) 


(4.7.13) 


Choosing P, _ -aP 2 one obtains, assuming P 2 > 1, a P 2 > 1: 


Wf.0)^(l + a)- 2 


(4.7.14) 


so that p may get close to 1 if a 
Four-direction point Gauss-Seidel 


is small. The remedy is to include more sweep directions, 
(consisting of four successive sweeps with four orderings: 
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£ 

p 

PD 

(3 

l 

0.25 

0.25 

0 

10' 1 

0.27 

0.25 

0 

10 -2 

0.45 

0.28 

105° 

10" 3 

0.71 

0.50 

105° 

10" 5 

0.77 

0.71 

105° 


Table 4.7.2: Fourier smoothing factors p,p D for the convection-diffusion equation discretized 
according to (4.5.11); symmetric point Gauss-Seidel smoothing. 

the forward and backward orderings of Figure 3.3.1, the forward vertical line ordering of Fig- 
ure 3.3.1, and this last ordering reversed) is robust for this test problem, as r ns y 

Table 4.7.3. 

As before, we have taken 0 = fcsr/12, t = 0, 1,2 23; Table 4.7.3 gives results only 

for a value of (3 for which the largest p and p D are obtained. Clearly, four- direct ion pom 
Gauss-Seidel is an excellent smoother for the convection-diffusion equation. It is found 

o and on change little when n is increased further. . 

Another useful smoother for this test problem is four-direction point Gauss-Seidel- Jacobi, 


£ 

p 

PD 

P 

i 

0.040 

0.040 

0° 

10" 1 

0.043 

0.042 

0° 

10“ 2 

0.069 

0.068 

0° 

10 -3 

0.16 

0.12 

0° 

10" 5 

0.20 

0.0015 

15° 


Table 4.7.3: Fourier smoothing factors p,p D for the convection-diffusion equation discretized 
according to (4.5.11); four-direction point Gauss-Seidel smoothing; n - 64. 


defined in Section 3.3. As an example, we give for discretization (4.5.11) the splitting for the 
forward step: 

1 0 1 

+ f[-c-M 2|c| 0] 


[M] = e 


-14 0 
0 

[N] = [M] - [A] 


(4.7.15) 
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The amplification factor is easily derived. Table 4.7.4 gives results, sampling 0 as before. The 
results are satisfactory, but there seems to be a degradation of smoothing performance in the 

vicinity of/? 0 and similarly near /? = ki r/2, k = 1,2,3). Finer samphng with intervals 
of l u gives the results of Table 4.7.5. 


This smoother is clearly usable, but it is found that damping improves performance still 
further. Numerical experiments show that w = 0.8 is a good vaJue; each step is damped 
separately. Results are given in Table 4.7.6. Clearly, this is an efficient and robust smoother 
tor the convection-diffusion equation, with u fixed atw = 0.8. Choosing u> = 1 gives a little 
improvement for e/h >0.1, but in practice a fixed value of w is to be preferred, of course. 


P Pd 0 


1 

0.130 

10- 1 

0.130 

10 -2 

0.127 

10“ 3 

0.247 

10~ 5 

0.509 

10 -8 

0.514 


0.130 0* 

0.130 45° 
0.127 45° 
0.242 15° 

0.494 15° 

0.499 15° 


Table 4.7.4: Fourier smoothing factors p,p D for the convection-diffusion equation discretized 
according to (4.5.11); four-direction point Gauss-Seidel-Jacobi smoothing; n = 64. 


£ 

n 

P 

0 

PD 

0 

10" 8 

64 

0.947 


0.562 

8° 

10" 8 

128 

0.949 

i° 

0.680 

5° 


Table 4.7.5: Fourier smoothing factors p,p D for the convection-diffusion equation discretized 
according to (4.5.1 1); four-direction point Gauss-Seidel-Jacobi smoothing. 


Line Gauss-Seidel 

For forward vertical line Gauss-Seidel we have 

\(6) = e' 9 ' [1 - Pi - |/\|)/ 2 ]/{ 4 + |Pj| + \P 2 \ -e**i[l + (p, + |/> 1 |)/ 2 ] 
-e ,(?2 [l + (|F 2 | - P 2 )/2] - e**[l + (P 2 + \P 2 |)/2]} (4.7.16) 
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£ 

Pi Pd 

P 

1.0 

0.214 

0 U 

10" 1 

0.214 

0° 

10 -2 

0.214 

45° 

10“ 3 

0.217 

45° 

10-5 

0.218 

45° 

10' 8 

0.218 

45° 


Table 4.7.6: Fourier smoothing factors, p,p D for the convection-diffusion equation discretized 
according to (4.5.11); four- direction point Gauss-Seidel-Jacobi smoothing with damping pa- 

rameter u = 0.8; n = 64. 

For Pi < 0, P 2 > 0 this gives |A(tt, 0)| = (1 - Pi)/(3 - Pi), which tends to 1 as |Pi| I - 
so that this smoother is not robust. Alternating line Gauss-Seidel is also not robust for this 
test problem. If P 2 < 0, Pi = aP 2 , a > 0 and |P 2 | > 1, \aP 2 \ > 1 then 


A(0,x/2)Sia/(l + a-i) 


(4.7.17) 


m tW lAtn tt/211 “ a/ffl + a) 2 + 1] 1/2 , which tends to 1 if a > 1. Symmetric (forward fol- 

vLu’ca, Ga^eide. for 

Table 4.7.7 presents some results. Again, n — 64 and (3 — kn/2, 1 ’ ’ ’ ’ 


£ 

p 

p 

PD 

P 

i 

0.20 

90° 

0.20 

90° 

10 _1 

0.20 

CO 

o 

o 

0.20 

90° 

10" 2 

0.20 

o 

o 

0.20 

90° 

10~ 3 

0.30 

0° 

0.26 

0° 

10' 5 

0.33 

0° 

0.0019 

75° 


Table 4.7.7: Fourier smoothing factors, p,p D for the convection-diffusion equation discretized 
according to (4.5.11); symmetric vertical line Gauss-Seidel smoothing; n - 


gives results only for the worst case in 0. 

We will not analyse these results further. Numerically we find that for /3 - 0 and 
p = (A(0, */2) = (1 + Pi)/(9 + 3 Pi) = 1/3. As £ | 0, p D depends on the value 

clear that we have a robust smoother. 


e <C 1 that 
of ne. It is 
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We may conclude that alternating symmetric line Gauss-Seidel is robust for both test 
problems, provided the mixed derivative is discretized according to (4.5.8). A disadvantage 
this smoother is that it does not lend itself to vectorized or parallel computing. 

The Jacob.-type methods discussed earlier and Gauss-Seidel with pattern orderings (white- 
black, zebra) are more favourable in this respect. Fourier smoothing analysis of Gauss-Seidel 
with pattern orderings is more involved, and is postponed to a later section. 

Exercise 4.7.1 Show that damped point Gauss-Seidel is not robust for the rotated anisotropic 
diffusion equation with c = 1, s = 0, with standard coarsening. 

Exercise 4.7.2 An Exercise 4.7.1, but for Gauss-Seidel-Jacobi method. 

4.8 Incomplete point LU smoothing 

For Fourier analysis it is necessary that [M] and [N] are constant, i.e. do not depend on the 
location in the grid. For the methods just discussed this is the case if [A] is constant. For 
incomplete factorization smoothing methods this is not, however, sufficient. Near the bound- 
ary o the domam[M] (and hence [N] = [M]-[A]) varies, usually tending rapidly to a con- 
stant stencil away from the boundaries. Nevertheless, useful predictions about the smoothing 
performance of incomplete factorization smoothing can be made by means of Fourier analysis. 
How this can be done is best illustrated by means of an example. 

Five-point ILU 

This incomplete factorization has been defined in Section 4.4, in standard matrix notation 
In Section 4.4 A was assumed to have a five-point stencil. With application to test problem 

( .5.6) in mind, A is assumed to have the seven-point given below. In stencil notation we 
have 


[A] = 


[D)i = 



'/ 9 


0 

=z 

Ci d q 

. w, = 

c Si 0 


a b 


a 


0 


“ ■* 

9 

“ 

0 Si 0 

. m = 

0 S{ q 


- 0 J 


0 


(4.8.1) 


with a - 0* vciaiuu. rui vi we nave tne recursion (3.4.12) 

s t = d- ag/Si_ e2 - cq/Si_ ei ( 4 . 8 . 2 ) 

where e 1= (1,0), e 2 = (0,1). Terms involving negative values of i a , a = 1 or 2, are to be 
repaced by zero. We will show the following Lemma. 
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Lemma 4.8.1. If 


a + c + d + q + g> 0, a,c,q,g< 0 , d> 0 


(4.8.3) 


then 


lim 

X 1 ,t2 — f OO 


Si = 6 = d/2 + [d 2 /4 - (a</ + cg)] 1/2 


(4.8.4) 


The proof is given in [141]. Note that (4.8.3) is satisfied if b = f - 0 and A is z K-mznx 
(Section 3.2). Obviously, is real, and 6 > d. The rate at which the limit is reached m (4.8.4) 
L studied in [141]. A sufficient number of mesh points awa, from the Wdar^ of the grid 
G we have approximately 6, = 6 , and replacing by 6 we obtain for [M] - [L][ ][ ]• 


[M] = 


eg / 6 g 

c d q 

a aq /6 


(4.8.5) 


and standard Fourier smoothing analysis can be applied. Equation (4.8.5) derived 
ily by nothing that in stencil notation (ABu) t = Vp k A{t,j)B{' + j, k)u i+i+k , so that 
A(i j)B(i+j,k) gives a contribution to C(i,j + fc), where C = AB\ by summing contribu- 
tions one obtains C{i,l). An explicit expression for C(i,l) is C(t,l) - ZjA(i,j)B(i+J, j), 
since one can write (Cu)i = S/SjA(f, j)B(i + jj - j)ui+i- 


Smoothing factor of five-point ILU . 

The modified version of incomplete factorization will be studied. As remark«l m [145J 
fication is better than damping, because if the error matrix N rs smah with a - - 0 it al o 
be small with a / 0. The optimum <r depends on the problem. A fixed a for aU pr 
be preferred. From the analysis and experiments in [145] and [147] and our own experiments 
it follows that a = 0.5 is a good choice for al] point- factorizations considered here and all 
problems. Results will be presented with <r = 0 and rr = 0.5. The modtfied versron of the 

recursion (3.4.12) for Sk is 

6 k = d- ag/h-l ~ cq/h-\ + <r{\oq/h-i ~ b \ + “ /I) (4.8.6) 

The limiting value 6 in the interior of the domain, far from the boundaries satisfies (4.8.6) 
with the subscripts omitted, and is easily determined numerically by the following recursion 


S k+i = d-{aq + cq)/ 6 h + a{\aq/ 6 k - b\ + \cg/ 6 k f |) 


(4.8.7) 


The amplification factor is given by 

A(0) = {(aq /6 - b)exp[i( 0 \ - #2)] + - /)exp[t(02 _ + a P)/ 

{aexp(-W 2 ) + aqexp[i{ 6 i - 82 )}/ & + cexp(-t9i) + d + crp 
+qexp(i 0 i) + cgexp[i( 6 2 - 0 i )}/6 + gexp(i 02 )} 
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where p = | aq/S - b\ + \cg/b - f |. 

Anisotropic diffusion equation 

For the (non-rotated 0 = 0°) anisotropic diffusion equation with discretization (4.5.6) we have 

9 ~ a = c - q = -£, d = 2 + 2e, b = f = 0, and we obtain: S = 1 + e + \2e(l + a)! 1 / 2 
and n ’ 


A(0) = [£cos(^ -e 2 )l6 + ael6]/ 

[l + £ + <7£/tf- £COS0! - COS0 2 + £COs(<9! -<9 2 )/<5] 

We wil] study a few special cases. For £ = 1 and u = 0we find in [141]: 

P = |A(tt/ 2, — 7r/3)| = (2V3 + VG- I) -1 ~ 0.2035 
The case £ = 1, a ^ 0 is analytically less tractable. For e < 1 we find in [141]: 

0 <a< 1/2: pa|A(jr, 0 ) = (l-<r)/(2tf-l + ff) 

1/2 < cr < 1 : p = |A(7 t/ 2 ) = <j/((t + 6) (4.8.11) 

0< ( r< 1/2: Pd — |A(7t,t)| = (1 — a)j{26 — 1 + <r + 8t 2 /2e) 

1/2 < <7 < 1 : p D S |A(7r/2,r)| = (<7 + r)(<r + £ + 6t 2 /2s) (4.8.12) 

where r^= 27 t/ti 2 . These analytical results are confirmed by Table 4.8.1. For example, for 
e = 10 , n 2 = 64 and a = 1/2 equation (4.8.12) gives p D £ 0.090, p £ 1/3. Table 4.8.1 

includes the worst case for 0 in the set {/l = Attt/ 12, A: = 0, 1,2, ...,23). Here we have another 
example showing that the influence of the type of the boundary conditions on smoothing 
analysis may be important. For the non-rotated anisotropic diffusion equation ((3 = 0° or 
0 = 90 ) we have a robust smoother both for a = 0 and cr = 1/2, provided the boundary 
conditions are of Dirichlet type at those parts of the boundary that are perpendicular to the 
direction of strong coupling. When 0 is arbitrary, five-point ILU is not a robust smoother 
with C7 _ 0 or <7 _ 1/2. We have not experimented with other values of <7, because, as it wil] 
turn out, there are other smoothers that are robust, with a fixed choice of cr, that does not 
depend on the problem. 


(4.8.9) 

(4.8.10) 
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e 

a 

P 

0 = 0°, 90° 

P 

0= 15° 

PD 

P = 0°, 90° 

PD 

P = 15° 

i 

0 

0.20 

0.20 

0.20 

0.20 

KT 1 

0 

0.48 

1.48 

0.46 

1.44 

1(T 2 

0 

0.77 

7.84 

0.58 

6.90 

10" 3 

0 

0.92 

13.0 

0.16 

10.8 

10~ 5 

0 

0.99 

13.9 

0.002 

11.5 

1 

0.5 

0.20 

0.20 

0.20 

0.20 

10" 1 

0.5 

0.26 

0.78* 

0.26 

0.78* 

10" 2 

0.5 

0.30 

1.06 

0.025 

1.01 

10" 3 

0.5 

0.32 

1.25 

0.089 

1.18 

IQ' 5 

0.5 

0.33 

1.27 

0.001 

1.20 


Table 4.8.1: Fourier smoothing factors, p,p D for the rotated anisotropic diffusion equation 
discretized according to (4.5.6); five-point ILU smoothing; n = 64. In the cases marked wi 

*, P = 45° 

Convection-diffusion equation , * 

Let us take ft = -aft, a > 0, ft > 0, where ft = ch/e, ft = sh/e Then we have for the 

convection- diffusion equation discretized according to (4.5.11): a - 1 r 2 , J » c 

i A - 4 -i- H + olPo o = -1 - aP 2 , 9 - -1- After some manipulation one finds that it 
a<l P 2 >1, ocp\ »’ l then A(tt/2, 0) - i as P 2 - oo. This is accordance with Table 4.8.2. 
The worst case obtained when 0 is varied according to 0 = kir /12, k - 0, 1 , 2 ,.. .,23 is ste . 
Clearly, five-point ILU is not robust for the convection-diffusion equation, at least for a - U 

and a = 0.5 

Seven-point ILU . A TT TT a 

Seven-point ILU tends to be more efficient and robust than five-point ILU. Assume 


[A] 


/ 9 

c d q 

a b 


(4.8.13) 


The seven-point incomplete factorization A = LD~'U-N discussed in Section 4.4 is defined 
in stencil notation as follows: 


[L)i = 


0 0 

7i 0 
a; Pi 


, [D t = 


0 0 
0 6i 0 
0 0 


, [U},= 


C* v 

0 Si Pi 
0 0 


(4.8.14) 
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£ 

p 

PD 

P 

P 

PD 

cr - 0 

<7 = 

0.5 

1 

0.20 

0.20 

0°~ 

0.20 

0.20 

10- 1 

0.21 

0.21 

0° 

0.20 

0.20 

10 -2 

0.24 

0.24 

120° 

0.24 

0.24 

10~ 3 

0.60 

0.60 

105° 

0.48 

0.48 

10“ s 

0.77 

0.71 

105° 

0.59 

0.58 


Table 4.8.2: Fourier smoothing factors p,p D for the convection-diffusion equation discretized 
according to (4.5.11); five-point ILU smoothing; n — 64 


We have, taking the limit i — oo, assuming the limit exists and writing lim,-^ a, = a etc., 

a = a, (3 = b — a M / S, j — c — a(/S , 

p = q - 0g/6, ( = / - jg/S, g = g (4.8.15) 

with 6 the appropriate root of 

6 = d - (ag + /?£ + 7 p,)6 + a(\/3p/8\ + |tC/^|) (4.8.16) 

Numerical evidence indicates that the limiting 8 resulting as t oo is the same as that for 
the following recursion, 


0o — by 70 — — d , f . iq — q^ Co = f 

0j + 1 = ^ “ a Pj i^j ? 7j+i = c ~ a (j/$j 
6 J+1 =d-(ag + /3 j+ i(j + 7 'j+lPj)/6j + a(\P J+l p J /8 : \ + |7 J+ i0/0|) 
Pj+i = <7 ~ fij+ig/Sj, Ci+i = / - ~fj+\g/6j 

For M we find M = LD~ l U = A + N, with 


[fV] = 


p 2 0 0 

o P3 

0 0 0 Pl 


Pi 

P3 


Pp/8, p 2 = 7C/O 

CT (|Pi| + IP 2 I 1 ) 


The amplification factor is given by 


(4.8.17) 


(4.8.18) 


= {Pa +Pi exp[i( 20 , - 9 2 )\ + p a exp[-*( 2 tf 1 - 0 2 )]}/ 

{a exp(— i# 2 ) + 6exp[t(0 2 - flj)] + Pl exp[i(20, - 0 2 )] + cexp(-^ 1 ) 
+d + p 3 + 9 exp(i0 1 ) + p 2 exp[— i(20j - 0 2 )] 

+/ exp[— 1(0, - 0 2 )] + g exp (»0 2 )} 
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Po"tr:»itttprc S aXr„1rI 1 e m discretized according to (4.5.9) we have s ym me.r y : 
M = 7 , C = 0, 9 = a, f = b, q = c, so that (4.8.19) becomes 

A(0) = [crp + pcos(20i - 0 2 }/ 

[acos62 + bcos{9 l -6 2 ) + ccos6i + d/2 + ap + pcos{2e 1 -62)] (4.8.20) 

with p = 

With rotation angle (3 = 90° and e < 1 we find in [141]: 


0 < c < 1/2 : |A(0, tt)| ~ 

1/2 > a > 1 : p cs |A(0,ir/2)| ~ 


_££_ 
5+trp 


(4.8.21) 


0 < a < 1/2 : p D s |A(<p, tt)[ ~ |(a - 1 + V)/f(2 + V/ 2e ) + a ~ (4.8.22) 

1/2 < a < 1 : Pd — |A((p,?r/2)| ~ |(<r + 2<p)/[5 2 (l + <P /&) + * - M\ 
with «p = 2 ,/n,. These results agree approximately with Table 4.8. * For example, for 
e = 10" 3 tii = 64 Equation (4.8.22) gives p D ~ 0.152 for a - 0, and p D - 0.10 
Table 4.8.3 includes the worst case for f3 in the set {(3 = kn/12, k = 0,1,2,..., }• qua '°" s 

(4 8 21) and (4.8.22) and Table 4.8.3 show that the boundary conditions may have an imp - 
K For rotation angle fi = 0 or 0 = 90°, seven-point 1LU is a good smoother for 


£ 

a 

P 

(3 = 0° 

P 

(3 = 90° 

P,P 


O 

o 
Q " 

PD 

(3 = 90° 

PD,fi 


1 

o 

0.13 

0.13 

0.13, 

any 

0.12 

0.12 

0.12, 

any 

IQ- 

0 

0.17 

0.27 

0.45, 

75° 

0.16 

0.27 

0.44, 

75° 

10 -2 

0 

0.17 

0.61 

1.35, 

75° 

0.11 

0.45 

1.26, 

75° 

10 -3 

o 

0.17 

0.84 

1.69, 

75° 

0.02 

0.16 

1.55, 

75° 

10 -5 

0 

0.17 

0.98 

1.74, 

75° 

io - 4 

0.002 

1.59, 

75° 

1 

0.5 

0.11 

0.11 

1.11, 

any 

0.11 

0.11 

0.11, 

any 

10 -1 

0.5 

0.089 

0.23 

o 

on 

o 

60° 

0.087 

0.23 

0.50, 

60° 

a /y A 

10 -2 

0.5 

0.091 

0.27 

0.77, 

O 

O 

CO 

0.075 

0.25 

0.77, 

60° 

10 -3 

0.5 

0.091 

0.31 

o 

bo 

60° 

0.029 

0.097 

0.82, 

60° 

10- 5 

0.5 

0.086 

0.33 

0.83, 

60° 

4 x 10~ 4 

10' 3 

0.82, 

60° 


Table 4.8.3: Fourier smoothing factors p,p D for the rotated anisotropic diffusion equation 
discretized according (4.5.6); seven-point ILU smoothing; n = 64 

the anisotropic diffusion equation. With a = 0.5 we have a robust smoother; finer sampling 
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of 0 and increasing n gives results indicating that p and p D are bounded away from 1. For 

some values of 0 this smoother is not, however, very effective. One might try other values of 

I" ;™ mish PD - A more efficient and robust ILU type smoother will be introduced shortly 
In [141] it is shown that for (3 = 45° and e < 1 


max 


<r + 1 


o + 1 


} 


(4.8.23) 


Hence, the optimal value of a for this case is a — 0.5, for which p ~ 1/3. 


Convection-diffusion equation 

Table 4.8.4 gives some results for the convection- diffusion equation. The worst case for 0 is 
t e set {0 — k-KjVl : k = 0, 1, 2, ...,23} is listed. It is found numerically that p <C 1 and 
PD < 1 when e < 1, except for 0 close to 0° or 180°, where p and p D are found to be 
larger than for other values of 0, which may spell trouble. We, therefore, do some analysis. 

umerically it is found that for £ < 1 and |s| < 1 we have p ~ |A(0,7r/2)|, both for a = 0 
and a = 1/2. We proceed to determine A(0 ,tt/ 2). Assume c < 0, s > 0; then (4.5.11) gives 

U 8 r« ~ • = ~ £ ’ d = 4£ ~ Ch + * h ' q = + hc ' f = °, 9 = ^^5 

(4.8.15) and (4.8.16) give, assuming e < 1 , \s\ < 1 and keeping only leading terms in £ and 
(£ + sh ) ch / 6 i 7 — -£■> P — ch, C — 0, 6 ~ (s — c)h, /»,-(£ + sh)c 2 /(s - c) 2 , p 2 = 0. 
Substitution in (4.8.19) and neglect of a few higher order terms results in 


A(0, tt/2) ~ 

where r = sh/e, so that 


(a- i)(r + 1) 


(r + 2)(1 — 2 tan /?) + <r(l + r) + t(l — 2 r tan 0) 


(4.8.24) 



<7 = 0 



<7 = 

0.5 


e 

p 

PD 

0 

p 

PD 

0 

i 

0.13 

0.12 

90° 

0.11 

0.11 

~w~ 

10-t 

0.13 

0.13 

90° 

0.12 

0.12 

0° 

10~ 2 

0.16 

0.16 

0° 

0.17 

0.17 

165° 

10~ 3 

0.44 

0.43 

165° 

0.37 

0.37 

165° 

10“ 5 

0.58 

0.54 

165° 

0.47 

0.47 

165° 


Table 4.8.4: Fourier smoothing factors p,p D for the convection-diffusion equation discretized 
according to (4.5.11); seven-point ILU smoothing; n = 64 

P 2 -(t+ 1) v + 1 )/{[(r + 2)(1 - 2 tan 0) + a( 1 + r)] 2 + (1 - 2r tan 0 ) 2 } (4.8.25) 
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hence, 


p 2 < (cr 2 + l)/(cr+ i) 2 (4.8.26) 

Choosing a = 1/2, (4.8.26) gives p < * 0.75, so that the smoother is robust With 

<7 = 0, inequality (4.8.26) does not keep p away from 1. Equation (4.8.25) gives, for <r 

limp=lA/5, lim p = (1 -4 tan /3 + 8tan 2 /J) _1/2 (4.8.27) 

T — *-0 * T -^°° 

This is confirmed by numerical experiments. With (T = 1/2 we have a robust smoother for 
the convection-diffusion equation. Alternating ILU, to be discussed shortly, may, however 
be more efficient. With <r = 0 , p < 1 except in a small neighbourhood of 0 - 0 and 
3 = 180°. Since in practice r remains finite, some smoothing effect remains. For examp e 
for 5 = 0.1 (0 ~ 174.3), h = 1/64 and £ = 10" 5 we have r ~ 156 and (4.8.27) gives p ~ 0.82. 
This explains why in practice seven-point ILU with <r = 0 is a satisfactory smoother for the 
convection-diffusion equation but (7 = 1/2 gives a better smoother. 


Nine-point ILU 

Assume 


Reasoning as before, we have 
0 0 0 
[L] 


[A] = 


7 

LO 


6 

a 


0 

P 


D = 


/ 

9 

p 





c 

d 

9 





z 

a 

b 





0 0 

0 ' 



' C 

V 

T 

0 S 

0 

y 

u = 

0 

6 

P 

0 0 

0 



0 

0 

0 

in [141], here interpreted 

as 


(4.8.28) 


(4.8.29) 


unknowns. The relevant solution of three equations may De oDiameu 
following recursion 

a 0 = a, Po = b, 7o = c, = d, p 0 = <7, Co = /. Vo = 9 

aj + 1 = a - zpj/Sj, (3j+ 1 = b - aj+iPj/Sj 
7j+l = C - (zr)j + a j+ iCj)/^ 
rij+i = {\flj+\Pj\ + l z Cjl + \Pj+iP\ 4" l7j+iCjl)/^j 
6j +i = d - {zp + a j+\Vj + Pj+i<i + 7j+i/ i j)/^i + an i+ 1 
Pj+i = 9 ~ ( a j+lP + Pj+iVj)/ a i+i 
Cj+i = / ~ lj+iVj/ s ii Vj + 1 = 5 - 7j+iP/*i+i 

For M we find M = LD~ l U = A + N, with 

0 0 0 

z( 0 mi 0 /?p 

0 0 0 0p 


(4.8.30) 


N =8 


(4.8.31) 
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(4.8.32) 


with n - Ij(I + \zQ\ + \f3p\ + |/?/x|. The amplification factor is given by 

\(0) = B(6)/{(B(0) + A(0)} 

where 


B (0) = { 7 C exp [i(e 2 - 2 <9, )J + zCexp (-2i9 x ) 
+flp exp (2*0!) + (3p exp [t'(20! - 0 2 ] + <r n }/6 

and 


A(6) - z exp [ i(6] + 6 2 )\ + aexp (~i0 2 ) + 6 exp[*( 0 i - 0 2 )] + C exp (-iO^) 

+d + q exp (* 0 , ) + / exp [*( 0 2 - 9 , )] + < 7 exp (i 0 2 ) + pexp [i( 0 , + 0 2 )] 

Anisotropic diffusion equation 

For the anisotropic diffusion equation discretized according to (4.5.6) the nine-point ILU fac- 
torization is identical to the seven-point ILU factorization. Table 4.8.5 gives results for the 
case that the mixed derivative is discretized according to (4.5.8). In this case seven-point ILU 
performs poorly. When the mixed derivative is absent (/? = 0 ° or /? = 90°) nine-point ILU 
is identical to seven-point ILU. Therefore Table 4.8.5 gives only the worst case for (3 in the 
set {/? _ k/ 27T, k = 0, 1, 2, ..., 23}. Clearly, the smoother is not robust for a = 0 . But also 
°r ° = 1 / 2 there are va lues of 0 for which this smoother is not very effective. For example 
with finer sampling of (3 around 75° one finds a local maximum of approximately p D = 0.73 

lui p — oO • 



(7=0 




a — 

1 

2 



£ 

p 

a 

PD 

a 

P 

(3 

PD 

(3 

1 

0.13 

any 

0.12 

any 

0.11 

any 

0.11 

any 

10- 1 

0.52 

75° 

0.50 

75° 

0.42 

75° 

0.42 

60° 

10~ 2 

1.51 

75° 

1.34 

75° 

0.63 

75° 

0.63 

75° 

io ~ 3 

1.87 

75° 

1.62 

75° 

0.68 

75° 

0.68 

75° 

10~ 5 

1.92 

75° 

1.66 

75° 

0.68 

75° 

0.68 

75° 


Table 4.8 5: Fourier smoothing factors p,p D for the rotated anisoptropic diffusion equation 
discretized according to (4.5.6), but the mixed derivative discretized according to (4 5 8 V 
nine-point ILU smoothing; n = 64 v ' 
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Alternating seven-point ILU , . , . f 

The amplification factor of the second part (corresponding to the second backward grid po 
ordering defined by (3.4.16)) of alternating seven-point ILU smoothing, with factors deno e 
by L,£),U, may be determined as follows. Let [A] be given by (4.8.13). The stencil repre- 
sentation of the incomplete factorization discussed in Section 3.4 is 


[L} = 


0 7 


0 0 


Co 

0 <5 a 

, [£>} = 

o 

1^5 

o 

, [U} = 

r) 6 0 

0 p 


o 

o 

1 


ft 0 


(4.8.33) 


In [141] it is shown that a, 0, ..., f, are given by (4.8.15) and (4.8.16), provided the following 
substitutions are made: 


a — *• <7, 6 — *- 6, c —> g, d —* d, q —* a, /—*}■> 9 

The iteration matrix is M = LD~ X U = A + N. According to [141], 


(4.8.34) 


[N] 


P 2 

0 0 0 
0 p3 0 
0 0 0 
Pi 


(4.8.35) 


with ft = fit/l, ft = 7</«. P3 = »(lftl + Iftl)- It fo1o» 5 th »‘ the am pUhcation factor X(*) 
of the second step of alternating seven-point ILU smoothing is given by 

X(6) = {p 3 + Pi exp [i{ 6 1 - 202 )] + Pi ex P l*( 2<? 2 - )]}/ 

{a exp (~i0 2 ) + 6 ex p [i(0i - 0 2 )] + cex P W + d + P 3 + 9 exp (*0) 

+/exp [—*(0i - 02)] + pexp (*0 2 ) + Pi exp [i(0i - 20 2 )] 

+p 2 exp [i(20 2 - 0i )} (4.8.36) 

The amplification factor of alternating seven-point ILU is given by A(0)A(0), with A(0) given 
by (4.8.19). 

Anisotropic diffusion equation 

Table 4 8.6 gives some results for the rotated anisotropic diffusion equation. The worst case 
for 3 in the set {(3 = for/12, k = 0, 1,2, ...,23} is included. We see that with cr - 0.5 we have 
a robust smoother for this test case. Similar results (not given here) are obtained when the 
mixed derivative is approximated by (4.5.8) with alternating nine-point ILU. 
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P 

PD 



e 


a 

3 = 0°, 90° 

3 = 0°, 90° 

Pi Pd 

3 

i 


0 

9 x 10" 3 

9 x 10 -3 

9 x 10 -3 

any 

10- 

-1 

0 

0.021 

0.021 

0.061 

30° 

10- 

-2 

0 

0.041 

0.024 

0.25 

45° 

10- 

-3 

0 

0.057 

3 x 10 -3 

0.61 

45° 

10- 

-5 

0 

0.064 

10 -6 

0.94 

45° 

1 


0.5 

4 x 10~ 3 

4 x 10 -3 

4 x 10 -3 

any 

10- 

1 

0.5 

0.014 

0.014 

0.028 

15° 

10- 

2 

0.5 

0.20 

0.012 

0.058 

45° 

10- 

3 

0.5 

0.026 

2 x 10' 3 

0.090 

45° 

10- 

5 

0.5 

0.028 

0 

0.11 

45° 


Table 4.8.6: Fourier smoothing factors p,p D for the rotated anisotropic diffusion equation 
discretized according to (4.5.6); alternating seven-point ILU smoothing; n = 64 

Convection-diffusion equation 

Symmetry considerations imply that the second step of alternating seven-point ILU smoothing 
has, for e < \,p ~ 1 for 3 around 90° and 270°. Here, however, the first step has p < 1. 
Hence, we expect the alternating smoother to be robust for the convection-diffusion equation 
This is confirmed by the results of Table 4.8.7. The worst case for 3 in the set {3 = kn/12 • 
k = 0,1, 2, ...,23} is listed. ' 

To sum up, alternating modified point ILU is robust and very efficient in all cases. The use 
of alternating ILU has been proposed in [97], 



<y = 0 


a = 0.5 


£ 

PiPD 

3 

Pi PD 

3 

1.0 

9 x 10 -3 

0° 

4 x 10~ 3 

0 ° 

10- 1 

9 x 10~ 3 

0° 

4 x 10~ 3 

0° 

10~ 2 

0.019 

105° 

7 x 10- 3 

0° 

10- 3 

0.063 

o 

LO 

o 

1— H 

0.027 

120° 

10~ 5 

0.086 

105° 

0.036 

105° 


Table 4.8.7: Fourier smoothing factors p, p D for the convection-diffusion equation discretized 
according to (4.5.11); alternating seven-point ILU smoothing; n = 64 
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Modification has been analyzed and tested in [65], [97], [83], [82], [145] and [147]. 


4.9 Incomplete block factorization smoothing 

For the smoothing analysis of incomplete block factorization we refer to [141]. We present 
some results. 

Anisotropic diffusion equation 

Tables 4.9.1 and 4.9.2 give results for the two discretizations (4.5.6) and (4.5.8) of the rotated 
anisotropic diffusion equation. The worst cases for f3 in the set {/ 3 = kn/\2, k = 0,1,. ..,23} 
are included. In cases where the elements of D do not settle down quickly to values indepen- 
dent of location as one moves away from the grid boundaries, so that in these cases Fourier 
smoothing analysis is not realistic. 


e 

II 

o 

o 

P 

(3 = 90° 

P,0 

II ^ 

o 

o 

PD 

(3 = 90° 

PD,P 

i 

0.058 

0.058 

0.058, any 

0.056 

0.056 

0.056, any 

10 _1 

0.108 

0.133 

0.133,90° 

0.102 

0.116 

0.116,90° 

10~ 2 

0.149 

0.176 

0.131,45° 

0.195 

0.078 

0.131,45° 

10' 3 

0.164* 

0.194 

0.157*, 45° 

0.025* 

0.005 

0.157*, 45° 

10“ 5 

0.141 

0.120 

0.166*, 45° 

0° 

0 

0.166*, 45° 


Table 4.9.1: Fourier smoothing factors p, pn for the rotated anisotropic diffusion equation 
discretized according to (4.5.6); IBLU smoothing; n = 64. The symbol * indicates that 
the coefficients do not become constant rapidly away from the boundaries; therefore the 
corresponding value is not realistic. 


Convection-diffusion equation 

Table 4.9.3 gives results for the convection-diffusion equation, sampling ft as before. 

If is clear that IBLU is an efficient smoother for all cases. This is confirmed by the multigrid 

results presented in [107]. 

4.10 Fourier analysis of white-black and zebra Gauss-Seidel smoothing 

The Fourier analysis of white-black and zebra Gauss-Seidel smoothing requires special treat- 
ment, because the Fourier modes V>(0) as defined in Section 4.3 are not invariant under these 
iteration methods. The Fourier analysis of these methods is discussed in detail in [112]. They 
use sinusoidal Fourier modes. The resulting analysis is applicable only to special cases of the 
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P P PD pD 

£ (3 = 0° /? = 90° (i = 0° p = 90° 

io - 1 
10 -2 
10~ 3 
10~ 5 


0.058 

0.058 

0.056 

0.056 

0.108 

0.133 

0.102 

0.116 

0.49 

0.176 

0.096 

0.078 

0.164* 

0.194 

0.025* 

5 x 10 

0.141* 

0.200 

0.000* 

0.000 


Table 4.9.2: Fourier smoothing factors p, p D for the rotated anisotropic diffusion equation dis 

n M The'svmh T* With miX ' d dCTiva,i ™ 

64. The symbol has the same meaning as in the preceding table. 
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p 

0 

PD 

a 

1.0 

0.058 

0° 

0.056 

~r 

10" 1 

0.061 

0° 

0.058 

0° 

10~ 2 

0.092 

0° 

0.090 

0° 

10 -3 

0.173 

0° 

0.121 

0° 

10 -5 

0.200 

0° 

io - 3 

15° 


to u . lhe — - 

Fourfer Ide?''” 5 defi " ed SeCti °” 4 ' 5 ' Theref ° re We WiU co " li "” e exponential 

The amplification matrix 

Speciabzing to two dimensions and assuming n, and n 2 to be even, we have 


i>j(0) = exp (ijO) 


with 

and 


) 

(4.10.1) 

• • 'Ha 1 

(4.10.2) 

— 771 a + 1 , . . . , 771 a + 1 } 

(4.10.3) 
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where m a = n a /2 — 1. Define 


€ e 5 = 0n[-7r/2,7T/2) 2 , 0 2 


= 0 1 - 


0 3 = 0 1 - 


0 

sijfn(02) 7r 


0 4 = 0 1 


( sign(6\)ir 
y sign{0\)-K 

( sign(0\)n ^ 
' 0 ) 


(4.10.4) 


, • - _1 f < 0- sionftl =1 t > 0. Note that 0 5 almost coincides with the set of 

smooth wavenumbers 0~ defined by (4.4.20). As we will see, Span {*(#'), V>(» 2 ). «« )} 

^rFouHt«p^“o"' of .n a rbit ra ry periodic 

grid function (4.3.7) can be written as 


u j = Y1 c 9 v»j(«) 


(4.10.5) 


with eg a vector of dimension 4. 

If the error before smoothing is cJV(«), «*" sm00tM ‘ lg il is 8 iven by (a(*)c#) T lW. 
with A(0) a 4 x 4 matrix, called the amplification matrix. 


TInT of* smooth wavenumbers 0, has been defined by (4.4.20). Comparison with ©i i as 
defined by (4.10.4) shows that ^(0 k ), k = 2,3,4 are rough Fourier modes, whereas 4>(6 ) 
is smooth, except when 0} = -»/2 or 0 4 = -n/2. The projection operator on the space 
spanned by the rough Fourier modes is, therefore, given by the following diagonal matrix 


/ 8 ( 0 ) 


\ 


Q(6) = 


(4.10.6) 


with «(*) = 1 if o, = — rr/2 and S 2 = -w/2, and «(») = 0 otherwise. Hence, a suitable 
definition of the Fourier smoothing factor is 

p = mAx{ X (Q(0)A(0)) :tfe0 5 } (4.10.7) 


with x the spectral radius. 

The influence of Dirichlet boundary conditions can be taken into account 
similar way as before. Wavenumbers of the type (0,0$) and (0f,O), s - 


heuristically in a 
1,3,4, are to be 
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disregated (note that 0 2 = 0 cannot occur), that is, the corresponding elements of c& are to 
be replaced by zero. This can be implemented by replacing QA by PQA with 


P(0) = 


( pi(0) \ 

i o 

0 p 3 {0) 

\ P4(0) } 


(4.10.8) 


where p x (9) = 0 if = 0 and/or 0 2 = 0, and pi(9) = 1 otherwise; p 3 (0) = 0 if 0\ - 0 (hence 
9\ = 0), and p 3 {9) = 1 otherwise; similarly, p${9) = 0 if 9 2 = 0 (hence 0 2 = 0), and p^{9) = 1 
otherwise. The definition of the smoothing factor in the case of Dirichlet boundary conditions 
can now be given as 

p D = max {x(P(0)Q(0)-<i(0)) : 9 £ ©g} (4.10.9) 


Analogous to (4.4.23) a mesh-size independent smoothing factor p is defined as 


p = sup{x(Q(#M(0)) : 9 £ 0 4 } (4.10.10) 


with 0 a = (-7t/2, 7t/2) 2 . 

White-black Gauss-Seidel 

Let A have the five-point stencil given by (4.8.1) with 6 = / = 0. The use of white- 
black Gauss-Seidel makes no sense for the seven-point stencil (4.8.1) or the nine-point stencil 
(4.8.28), since the unknowns in points of the same colour cannot be updated independently. 
For these stencil multi-coloured Gass-Seidel can be used, but we will not go into this. 

Define grid points (j\,h) with j\ + j 2 even to be white and the remainder black. We 
will study white-black Gauss-Seide with damping. Let e° be the initial error , e 1 / 3 the error 
after the white step, e 2 / 3 the error after the black step, and e 1 the error after damping with 
parameter u>. Then we have 

e) /3 = -«_ e2 + ce?_ e| +qe% ej +ge° J+e2 )/d, j, + j 2 even (4.10.H) 

e] /3 = £?, ji + h odd. 

The relation between e 2 / 3 and e 1 / 3 is obtained from (4.10.11) by interchanging even and odd. 
The final error e 1 is given by 

e) = ue)' 3 + {\ - u)e) (4.10.12) 

Let the Fourier representation of e a , a = 0, 1/3, 2/3, 1 be given by 

e] = £ cfW) . 

see. 
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(4.10.13) 


If e° 3 = ipj{0 3 ), s = 1,2,3 or 4, then 


s 1 / 3 = n(O s )il)j(0 a ), ji + j 2 even 
e) /3 = i’ J (0 s ), h+j 2 odd 

with fi(0) = — [aexp (—i0 2 ) + cexp (~i0\) + < 7 exp (i6\) + gex p (i0 2 )/d. Hence 


so that 


e) 13 = \(h( 0 s ) + l)exp (ijP) + ±(n(0 s ) - 1 ) 
x exp [iji(0? - 7r)J exp [ij 2 (6 2 - 7r)] 



1 

2 


( 1 + Mi 
Mi ~ 1 
0 
0 


-1 - Hi 0 

1 - Hi 0 

0 1 + Hi 

0 n 2 -\ 


0 

0 

1 Hi 
1 - Hi 


\ 


c °e 


} 


0 € 0 5 


(4.10.14) 


(4.10.15) 


where ^ = /i(< 9), fi 2 = (aexp (- i0 2 ) - cexp (- 1 O 1 ) - qex p (i0 1 ) + ^exp ( id 2 ))/d . If the 
black step is treated in a similar way one finds, combining the two steps and incorporating 
the damping step, 

c l = WA(9) + (1 — w)/}c® (4.10.16) 

with 


m = j 


( MiO 4- Hi ) 
MiO - Hi) 
0 
0 


Hence 


— Mi(l + Mi) 0 

Mi(Mi - 1) 0 

0 M 2 (l + Hi) 

0 M 2 (1 M 2 ) 


0 \ 
0 

— M2(l + Hi) 
Hi(Hi - 1) 


(4.10.17) 


P(9)Q(0)A(Q) 



^ Mi^MiO + Mi) 

-pi^Mi(l + Mi) 

0 

0 \ 


1 

Mi(l - Mi) 

Mi (Mi “ 1) 

0 

0 


“ 2 

0 

0 

P3Hi(1 + Hi) 

-P 3 M 2(1 + M 2 ) 



0 

0 

P 4 M 2(1 - M 2 ) 

PiHiiHi - 1) / 

(4.10.18) 

The eigenvalues of PQA are 






Ai(0) — 0, \ 2 (^) — ^MiiMi — 

1 + Mi^(l + Mi)}> 



A 3 (0) = O, A 4 (0) = ^Hib>3 -Pa + Hi(P3 + p 4 )} (4.10.19) 


68 



and the two types of Fourier smoothing factor are found to be 

p,p D = max {|wA 2 (0) + 1 - w|, \uj\ 4 (8) + 1 - w | : 8 £ 0 5 } (4.10.20) 

where p\ = p 3 = p 4 = 1 in (4.10.19) gives p, and choosing pi, p 3 , p 4 as defined after equation 
(4.10.8) gives pp. 


With u> = 1 we have p = p D = 1/4 for Laplace’s equation [112], This is better than 
lexicographic Gauss-Seidel, for which p = 1/2 (Section 4.7). Furthermore, obviously, white- 
black Gauss-Seidel lends itself very well for vectorized and parallel computing. This fact, 
combined with the good smoothing properties for the Laplace equation, has led to some of 
the fastest Poisson solvers in existence, based on multigrid with white-black smoothing [131 
[113]. 


Convection-diffusion equation 

With p = 0 equation ( 4 . 5 . 11 ) gives a = -e, c = -e - h, d = 4e + h, q = -e, g = -e, so 
that /xi,2(0, — jt/2) = (2 + P)/( 4 + P), with P = h/e the mesh Peclet number. Hence, with 
Pi = P3 = Pa = 1 we have A 2 , 4 (o, - jt/2) = (2 + P) 2 /( 4 + P) 2 , so that p -+ 1 as P -+ oo for all 
w, and the same is true for pp. Hence white-black Gauss-seidel is not a good smoother for 
this test problem. 


Smoothing factor of zebra Gauss-Seidel 

Let A have the following nine-point stencil: 


[A} = 


f 9 P 
c d q 
z a b 


(4.10.21) 


Let us consider horizontal zebra smoothing with damping. Define grid points (jfi,jf 2 ) with j 2 
even to be white and the remainder to be black. Let e° be the initial error, e 1 / 3 the error 
after the ‘white’ step, e 2 / 3 the error after the ‘black’ step, and e 1 the error after damping 
with parameter u Then we have 


c£ )-e l + de ) /3 + qe l / + 3 ei 

= -e 2 + a£< j-e 2 + ^j+e, -e 2 + f £ °j-e 1+ e 2 + P £ j+e 2 + P^j+e, +e 2 ), 



j 2 even 

ji odd (4.10.22) 

(4.10.23) 


where ej = ( 1 , 0) and e 2 = (0, 1 ). 
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The relation between e 2 / 3 and e 1/3 is obtained from (4.10.23) by interchanging even and 
odd, and the final error e 1 is given by (4.10.12). 

It turns out that zebra iteration leaves certain two-dimensional subspaces invariant in Fourier 
space. In order to facilitate the analysis of alternating zebra, for which the invariant subspaces 
are the same as for white-black, we continue the use of the four-dimensional subspaces V>(0) 
introduced earlier. 

In [141] it is shown that the eigenvalues of P(8)Q(9)A(8) are 

Aj(0) = 0, A 2 (0) = ^Pi^Pi(l "b Pi) — 2 P3P\{y — ^3(0) - 0 (4 10 24) 

A 4 (0) = |p2(l + P 2 ) + \piPi{P-2 — 1) 


with 


6 ) = - {z exp( — z(0i + 0 2 )] + a exp(-t0 2 ) + & exp[z'(0i - #2)] 

+ / exp[i(0 2 - 01 )] + 9 exp(0 2 ) + p exp[i(0i + 0 2 )}/ 

[c exp(— ?0i) + d + q exp(i0i)] 

and p 2 = Pi(0i — 7T, 02 - 7r). 

The two types of Fourier smoothing factor are given by (4.10.20), taking A 2 ,A 4 from 
(4.10.24). 

Anisotropic diffusion equation 

For £ = 1 (Laplace’s equation), u> = 1 (no damping) and p\ = P3 — P4 = 1 (periodic boundary 
conditions) wehavepi(0) = cos 02/(2 — cos 8\) and p 2 (0) = — cos 02/(2 + cos 0i) . One finds 
max {|A 2 (0)| : 0 € ©5} = |A 2 (7r/2,0)| = \ and max {|A 4 (0)| : 0 6 ©5} = |A 4 (7r/2, 7r/2)| = ^,so 
that the smoothing factor is p = p = ^. 

For e < 1 and the rotation angle (3 = 0 we have strong coupling in the vertical direction, so 
that horizontal zebra smoothing is not expected to work. We have p 2 (0) = — cos 0 2 /(l + £ + 
£ cos 0i), so that |A 4 (ir/2, 0)| = (1 + £) -2 , hence lim p D > 1. Furthermore, with <p = 2w /n, 

we have |A 4 (tt/ 2), <p)| = cos 2 <p/(l + £) 2 , so that limpo > 1 - 0(h 2 ). Damping does not help 

here. We conclude that horizontal zebra is not robust for the anisotropic diffusion equation, 
and the same is true for vertical zebra, of course. 

Convection-diffusion equation 

With convection angle (3 = 7r/2 in (4.5.11) we have 

/z 2 (0) = [(1 + P) exp ( -i0 2 ) + exp {i0 2 )\/(4 + P + 2 cos 6 X ) , 
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where P = h/e is the mesh Peclet number. With p 4 = 1 (periodic boundary conditions) we 
have A 4 - p 2 , so that A 4 (tt/2, 0) = (2+ P) 2 /(4 + P) 2 , and we see that u;A 4 (tt/2 , 0) + 1 - w »l 

for P 1, so that p ~ 1 for P ^ 1 for aD u>. Hence, zebra smoothing is not suitable for the 
convection-diffusion equation at large mesh Peclet number. 


Smoothing factor of alternating zebra Gauss-Seidel 

As we saw, horizontal zebra smoothing does not work when there is strong coupling (large 
diffusion coefficient or strong convection) in the vertical direction. This suggests the use of 
alternating zebra: horizontal and vertical zebra combined. Following the suggestion in [1121, 
we will arran ge alternating zebra in the following ‘symmetric’ way: in vertical zebra we do 
first the black’ step and then the ‘white’ step, because this gives slightly better smoothing 

factors, and leads to identical results for 0 = 0° and 0 = 90°. The 4 x 4 amplification matrix 
of vertical zebra is found to be 


where 


A v {8) 


1 

2 


V\{vi + 1 ) 
0 
0 


0 

v 2 {v 2 + 1) 
Va(va - 1 ) 
0 


0 ^1(^1 + 1 ) ^ 

^2(^2 +1) 0 

^2(^2 - 1) 0 

0 u x (v x - 1) ) 


(4.10.25) 


= ~( z ex P H(0 1 + # 2 )] + bexp [i(9 x - 0 2 )] + cex p (-i9 x ) 

+9 exp (*0j) + / exp [i(0 2 - 0i)] + pexp [i{6 x + 0 2 )}} / 

[a exp (~i9 2 ) + d + g exp (i0 2 )] 

and v 2 (9) = u 1 (9 1 -n, 9 2 -tc). We will consider two types of damping: damping the horizontal 
and vertical steps separately (to be referred to as double damping) and damping only after 
the two steps have been completed. Double damping results in an amplification matrix given 

A = PQ[( 1 - u d )I + u> d A v ][(l - u d )I + uj d Ah] (4.10.26) 

where A h is given in [141]. In the case of single damping, put u; d = 1 in (4.10.26) and replace 

A : (1 — u s )I + u s A (4.10.27) 

The eigenvalues of the 4 x 4 matrix A are easily determined numerically. 

Anisotropic diffusion equation 

Tables 4.10.1 and 4.10.2 give results for the smoothing factors p,p D for the rotated anisotropic 
iffusion equation. The worst cases for the rotation angle 0 in the set {0 = kir/\2, k = 
0, 1,2,. ..,23} are included. For the results of Table 4.10.1 no damping was used. Introduction 
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of damping (u d /lorw^l) gives no improvement. However, as shown by Table 4.10.2, 
the mixed derivative is discretized according to (4.5.8) good results are obtained. For cases 
with e = 1 or P = 0° or P = 90° the two discretizations are identical of course, so for these 
cases without damping Table 4.10.1 applies. For Table 4.10.2 P has been sampled with an 
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Table 4.10.1: Fourier smoothing factors p,p D for the rotated anisotropic diffusion equation 
discretized according to (4.5.6); alternating zebra smoothing; n = 64 
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0.229 

30° 
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0.503 

8° 
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0.653 
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0.537 
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0.668 

8° 
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0.538 

4° 

0.300 

0.668 
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Table 4.10.2: Fourier smoothing factors p,po for the rotated anisotropic diffusion equation 
discretized according to(4.5.6) but with the mixed approximated by (4.5.8); alternating zebra 
smoothing with single damping; n = 64 

interval of 2°. Symmetry means that only p 6 [0°,45°] needs to be considered. Results with 
single damping (u> s = 0.7) are included. Clearly, damping is not needed in this case and 
even somewhat disadvantageous. As will be seen shortly, this method, however, works for the 
convection diffusion test problem only if damping is applied. Numerical experiments show 
that a fixed value of = 0.7 is suitable, and that there is not much difference between single 
damping and double damping. We present results only for single damping. 
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Convection-diffusion equation 


“ e Ta^^“ n ,r™t d With in '\ rValS ° f2 “ ; the worsl cases are presented. The 

ow at alternating zebra without damping is a reasonable smoother 
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lt' e dhgto : (4Tnvr° 0, p” S fa n‘ 0rS " fOT the c°"v«ttio„-diffusio n equation discretized 
according to (4.5.11), alternating zebra smoothing with single damping; n = 64 

for the convection-diffusion equation. If the mesh Peclet numbers h cos 0/e or hsm 3/e 
ecomes arge ( > 100, say), p approaches 1, but p D remains reasonable. 

xe amping parameter u s = 0.7 gives good results also for p. The value u - 07 
chosen after some experimentation. 3 ~ U>7 


was 


We see that with u> s 0.7 alternating zebra is robust and reasonably efficient for both 

convection-diffusion and the rotated anisotropic diffusion equation, provided the mixed 
derivative is discretized according to (4.5.8). P ne mixed 

4.11 Multistage smoothing methods 

m “ lt ,“ ta f *r« thin * methods are «!«> of the basic iterative method type 

it uilrtoked " 6 Wi “ be “ PlaiMd) ' b “‘ ia — 'figHd literature rtey 

usually looked upon as techniques to solve systems of ordinary differential equations 

dXemLTeqltrf 1 diSCretoti °" “ f "■*«“ hyperbolic or almost hyperbolic partial’ 

The convection-diffusion test problem (4.5.4) is of this type, but (4.5.3) is not We will 
herefore consider the application of multistage smoothing to (4.5.4) only. Multistage meth’ 
ods have been introduced in [74) for the solution of the Euler equations o^^dinamL and 
as smoothing methods in a multigrid approach in [71). For the simple scalar test problem 
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(4 5 4) multistage smoothing is less efficient than the better ones of the smoothing methods 
discussed before" The simple test problem (4.5.4), however, lends itself well for explammg 
the basic principles of multistage smoothing, which is the purpose of this section. 

Artificial time- derivative , , 

The basic idea of multistage smoothing is to add a time- derivative to the equation to be 
solved, and to use a time-stepping method to damp the short wavelength components of the 
error The time-stepping method is of multistage (Runge-Kutta) type. Damping o 
waves occurs only if the discretization is dissipative, which implies that for hyperbohcor 
almost hyperbolic problems some form of upwind discretization must be used, or an art ficial 
dissipation term must be added. Such measures are required anyway to obtain good soluti . 
The test problem (4.5.4) is replaced by 


r\ 

— — e{u 11 + ^, 22 ) + cu y 1 + — / 

dt 


(4.11.1) 


Spatial discretization according to (4.5.10) or (4.5.11) gives a system of ordinary differential 
equations denoted by 

^ = -h~>Au + f (4.11.2) 

dt 

where A is the operator defined in (4.5.10) or (4.5.11); « is the vector of grid function values. 


The time-derivative in (4.11.2) is an artefact; the purpose is to solve Au -hf. Hence, e 
temporal accuracy of the discretization is irrelevant. Denoting the time-level by a superscrip 
n and stage number k by a superscript (fc), a p-stage (Runge-Kutta) discretization of (4.1 . ) 

is given by 


u {0) = u n 

u (k) _ u (°) - Ck „h- l Au( k - l) + c k Atf, k = l,2,...,p 

u n + i = (4.11.3) 

with c p = 1. Here v = At/h is the so-called Courant-Frederichs-Lewy (CFL) number. Elimi- 
nating u( k \ this can be rewritten as 

u n+l = P p {-vh~' A)u n + Qp-ii-vh,- 1 A)f (4.11.4) 


with the amplification polynomial P v a polynomial of degree P, defined by 

Pp(z) = 1 + z(l + Cp_iz(l -1- Cp_ 2 -z(...(l + Ciz)...) (4.11.5) 

and Q p -i is polynomial of degree p - 1 which plays no role in further discussion. 
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Semi-iterative methods 

Obviously, equation (4.1 1.6) can be interpreted as an iterative method for solving h~ 2 Au = f 
of the type introduced in Section 4.1 with iteration matrix 

S=P p {-vh~ x A) (4.11.6) 

Such methods, for which the iteration matrix is a polynomial in the matrix of the system to 
be solved, are called semi-iterative methods. See [129] for the theory of such methods. For 
P = 1 (one-stage method) we have 


S = I — vh l A (4.11.7) 

which is in fact the damped Jacobi method (Section 4.3) with diagonal scaling ( diag (A) = /), 
also known as the one-stage Richardson method. As a solution method for differentia] equa- ? 
tions this is known as the forward Euler method. Following the trend in the multigrid liter- 
ature, we will analyse method (4.11.3) as a multistage method for differentia] equations, but 
the analysis could be couched in the language of linear algebra just as well. 

The amplification factor 

The time step A t is restricted by stability. In order to assess this stability restriction and the 
smoothing behaviour of (4.11.4), the Fourier series (4.3.7) is substituted for u. It suffices to 
consider only one component u = 0 <= 0. We have vh~ l Axb(6) = uhr' utOWOY With 

A defined by (4.5.11) one finds 

p(0) = 4e + h{\c\ + |s|) - (2s + h|c|) cos 

— (2f + /i|s|) cos #2 + ihc sin + ihs sin 62 

and 

u n+1 = g(0)u n 

with the amplification factor g(6) given by 

9 ( 0 ) = P p (-ufi(e)/h) (4.11.10) 

The smoothing factor 

The smoothing factor is defined as before: 

P = max {|$(0)| :6e O r } (4.11.11) 

in the case of periodic boundary conditions, and 

p D = max{|p(0)| : 9 6 0?} (4.11.12) 


(4.11.8) 

(4.11.9) 
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for Dirichlet boundary conditions. 


Stability condition 

Stability requires that 

|flf(*)| < 1, V0 € 0 

The stability domain D of the multistage method is defined as 

D = {zeC : | P p (z) < 1} 


(4.11.13) 

(4.11.14) 


Stability requires that v is chosen such that 2 = -vfi(0)/h € D, V0 G 0- If P < 1 ut 
(4.11.13) is not satisfied, rough modes are damped but smooth modes are amplified, so that 

the multistage method is unsuitable. 


Local time-stepping , , , 

When the coefficients c and s in the convection-diffusion equation (4.11.1) are replaced by 

general variable coefficients v, and v 2 (in fluid mechanics applications v u v 2 are fluid velocity 
components), an appropriate definition of the CFL number is 

v = vAt/h, v = |ui| + |u 2 | (4.11.15) 

Hence, if A t is the same in every spatial grid point, as would be required for temporal ac- 
curacy will be variable if v is not constant. For smoothing purposes it is better to fix v 
at some favourable value, so that At will be different in different grid points and on different 
grids in multigrid applications. This is called local time-stepping . 

Optimization of the coefficients , 

The stability restriction on the CFL number y and the smoothing factor p depend on the 
coefficients c k . In the classical Runge-Kutta methods for solving ordinary differential equa- 
tions these are chosen to optimize stability and accuracy. For analyses see for example [115], 
[106]. For smoothing c k is chosen not to enhance accuracy but smoothing; smoothing is also 
influenced by v. The optimum values of v and c k are problem dependent. Some analysis 
of the optimization problem involved may be found in [127]. In general, this optimization 
problem can only be solved numerically. 


We proceed with a few examples. 

A four-stage method . 

Based upon an analysis of Catalano and Deconinck (prive-communication), in which optimal 
coefficients c k and CFL number v are sought for the upwind discretization (4.5.11) of (4.11.1) 
with e = 0, we choose 

ci = 0.07, C 2 = 0.19, c 3 = 0.42, u = 2.0 (4.11.16) 
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€ 

o 

o 

II 

0 = 15° 

~0 = 30° 

0 = 45° 

0 

10~ 5 

1.00 

0.997 

0.593 

0.591 

0.477 

0.482 

0.581 

0.587 


Table 4.11.1: Smoothing factor p for (4.11.1) discretized according to (4.5.11); four-stage 
method; n = 64 


Table 4.11.1 gives some results. 

Jt ‘ S found that P° differs ver y Uttle from p. It is not necessary to choose 0 outside 
[0 ,45 ], since the results are symmetric in 0. For £ £ 10~ 3 the method becomes unstable 
for certain values of 0. Hence, for problems in which the mesh Peclet number varies widely 
m the domain it would seem necessary to adopt c k and v to the local stencil. With e = 0 all 
multistage smoothers have p = 1 for grid-aligned flow (0 = 0° or 90°) : waves perpendicular 
to the flow are not damped. 

A five-stage method 

The following method has been proposed in [73] for a central discretization of the Euler 
equations of gas dynamics: 


ci - 1/4, c 2 - 1/6, c 3 = 3/8, c 4 = 1/2 (4.11.17) 

The method has also been applied to the compressible Navier-Stokes equations in [75], We 
will apply this method to test problem (4.11.1) with the central discretization (4.5.10) Since 
p(0) = ih(c sin + 5 sin 0 2 ) we have /x(0,tt) = 0, hence | ff (0,ir)| = 1, so that we have no 
smoother. An artificial dissipation term is therefore added to (4.11.2), which becomes 


with 


du 

dt 


— —h 2 Au — h 1 Bu -f- f 


[B] = x 


1 

-4 

1 -4 12 -4 1 

-4 
1 


(4.11.18) 


(4.11.19) 


where x is a parameter. 

We have B^Q) - r/(0)V>(0) with 
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“o 0- 

1 5° 

30° ~ 

“45^ 

p 

0.70 

0.77 

0.82 

0.82 


Table 4.11.2: Smoothing factor p for (4.11.1) discretized according to (4.5.10), five stage 
method; n = 64 


r,(0) = 4x[(l - cos ffj) 2 + (1 - cos 0 2 ) 2 } (4.11.20) 


For reasons of efficiency the artificial dissipation term is updated in [73] only in the first two 
stages. This gives the following five-stage method. 

U W = n(°) - c k v(h~ l A + B)ul k -'\ k = 1,2 (4.11.21) 

U W = — Ckv(h~^ k = 3,4,5 

The amplification polynomial now depends on two arguments zr,* 2 defined by z x = vh~'p{6), 
z 2 - vi ](0), and is given by the following algorithm: 


P 1 = 1 - C\{Z\ + Zg), P2 = 1 - c l( z l + z 2 )P\ 

P 3 = 1 - C 3 Z t P 2 - C 3 Z\ P 2 - C 3 Z 2 Pl, P4 = 1 - C4*l P3 - C 4 Z 2 Pl 
p 5 (z\,z 2 ) = 1 - Z\P 4 - Z 2 P\ 


(4.11.22) 


In one dimension Jameson and Baker [73] advocate «/ = 3 and X = 0.04; for stability </ should 

not be much larger than 3. In two dimensions max {vh ~ v( c + s ) - f ' “ S '“ 

Vy /2 = 3 gives 1/ ~ 2.1. With i/ = 2.1 and X = 0.04 we obtain the results of Table 4.11 2, for 

both e = 0 and e = 10" 5 . Again, p D ~ p. This method allows only £ < 1; for example, for 

£ - 10 -3 and ft = 45° we find p = 0.96. 
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Final remarks 

Advantages of multistage smoothing are excellent vectorization and parallelization poten- 
t.al, and easy generalization to systems of differential equations. Multistage methods are in 
wi esprea use for hyperbolic and almost hyperbolic systems in computational fluid dynam- 
ics They are not, however, robust, because, like all point-wise smoothing methods, they do 
not work when the unknowns are strongly coupled in one direction due to high mesh aspect 
ratios Also their smoothing factors are not small. Various strategems have been proposed in 
the Literature to improve multistage smoothing, such as residual averaging, including implicit 
stages, and local adaptation of c*, but we will not discuss this here; see [73], [75] and [127], 

4.12 Concluding remarks 

In this chapter Fourier smoothing analysis has been explained, and efficiency and robustness 
of a great number of smoothing methods has been investigated by determining the smoothing 
tactors p and p D for the two-dimensional test problems (4.5.3) and (4.5.4). The following 

met o s work for both problems, assuming the mixed derivative in (4.5.3) is suitably dis- 
cretized, either with (4.5.6) or (4.5.8): * 

(i) Damped alternating Jacobi; 

(ii) Alternating symmetric line Gauss-Seidel; 

(iii) Alternating modified incomplete point factorization; 

(iv) Incomplete block factorization; 

(v) Alternating damped zebra Gauss-Seidel. 

Where damping is needed the damping parameter can be fixed, independent of the problem. 

is important to take the type of boundary condition into account. The heuristic way in 
w ich this has been done within the framework of Fourier smoothing analysis correlates well 
with multigrid convergence results obtained in practice. 

Generalization of incomplete factorization to systems of differential equations and to non- 
linear equations is less straightforward than for the other methods. Application to the incom- 
pressible Navier-Stokes equations has, however, been worked out in [1441, [1461 [1481 ri 50 l 
and [149], and is discussed in [141]. ’ ’ 

Of course in three dimensions robust and efficient smoothers are more elusive than in two 
dimensions. Incomplete block factorization, the most powerful smoother in two dimensions 
is not robust in three dimensions [81]. Robust three-dimensional smoothers can be found 
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among methods that solve accurately in planes (plane Gauss-Seidel) [114] For a rooea»M 
multigrid approach to a complicated three-dimensional problem using ILU type smoothing, 

see [124], [122], [125], [123]. 

5 Prolongation, restriction and coarse grid approximation 

5.1 Introduction 

In this chapter the transfer operations between fine and coarse grids are discussed. 

ThTdom^n n in which the partial differential equation is to be solved is assumed to be the 
d-dimensional unit cube. In the case of vertex-centered discretization, the computational grid 

is defined by 

G - {x € IR d : x = jh, j = h = (huh 2 ,...,h d ), 

j a = 0,1,2 ,...,n a , h a = l/n a , a = l,2,...,d} ( 5,L1 ) 

In the case of cell-centered discretization, G is defined by 

G = {x £ R d : x = (j - s)h, j = (ji s = (l,l,...,l)/2, 

h = (h U h 2 ,";h d ), j a = 1,2, ..., n a , ha = 1 /»«, « = 1 » 2 — d > ( 5 - 1 ’ 2 ) 

These grids, on which the given problem is to be solved, are called fine grids Without danger 
of confusion, we will also consider G to be the set of d-tuples j occunng m (5.1.1) or (5.1.2). 

I^ThiTchapter it suffices to consider only one coarse grid. From the vertex- centered grid 
(5 1 1) a coarse grid is derived by vertex-centered, coarsening , and from the cell- centered gn 
(5 1 2 a coarse grid is derived by cell-centered coarsening. Coarse grid quantities will be 
identified by an overbar. Vertex- centered coarsening consists of deleting every other vertex in 
each direction. Cell-centered coarsening consists of taking unions of fine grid cells to obtain 
coarse grid cells. Figures 5.1.1 and 5.1.2 give an illustration. It is assumed that n a in (5.1. ) 

and (5.1.2) is even. 

Denote spaces of grid function by U : 


U = {w : G -*■ IR), U = {u:G->lR} 
The transfer operators are denoted by P and R: 

P : 0 —> U, R: U —* U 


(5.1.3) 


(5.1.4) 
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12 3 4 

I 1 1 I . 


*l-° 


*1 = 1 


* 1 "° 


*1 = 1 


Figure 5.1.1: Vertex-centered and cell-centered coarsening in one dimension. (• grid points) 
P is called prolongation , and R restriction. 

[141^ ere ° nly vertex ' centere( ^ coars ening will be discussed. For cell-centered coarsening, see 


5.2 Stencil notation 

In order to obtain a concise description of the transfer operators, stencil notation will be used. 

Stencil notation for operators of type U — ► U 

Let A : U — U be a Unear operator. Then, using stencil notation, Au can be denoted by 

(Au), = A(i,j)ui+j, ieG (5.2.1) 

with 7Z = {0,±1,±2, ...}. The subscript i = (t lt i 2 , ..., i d ) identifies a point in the computa- 
tional grid in the usual way; cf. Figure 5.1.2 for the case d = 2. 

The set Sa defined by 


S A ~ 0 € 2Z d : 3i 6 G with A(i,j) / 0} (5.2.2) 

is caUed the structure of A. The set of values A(i,j) with j € S A is caUed the stencil of A 
at grid point i. Often the word ’stencil’ refers more specifically to an array of values denoted 
by [AJj in which the values of A(i,j) are given; for example, in two dimensions, 


[A]i 


A{i, — e\ + e?) A(i,e 2 ) 
A(i, — €j) A(t, 0) 

Mh -c 2 ) 


A(i,e 1 ) 
A{i, e, - e 2 ) 


(5.2.3) 


where e, = (1,0) and e 2 = (0,1). 
[141]. 


For the representation of three-dimensional stencils, see 
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Vertex-centred 

Figure 5.1.2: Vertex- centered and cell-centered 


Cell-centred 

coarsening in two dimensions. 


(• grid points.) 


Example 5.2.1 One-dimensional discrete Laplacian: 

[A], = /r 2 [-i 2 -l] 


(5.2.4) 


Stencil notation for restriction operators 

Let R : U -* U be a restriction operator. Then, using stencil notation, Ru can be represente 

Y (Ru)i = R(i,j) u u+ji ie( 3 (5.2.5) 

jeZ d 


Example 5.2.2 Consider vertex-centered grids G,G for d - 1 as defined by (5.1.1) and as 
depicted in Figure 5.1.1. Let R be defined by 


Rui = WiU2i-\ + 2 U2i e * W2 *+ 1 ’ * — n 


(5.2.6) 


82 







with w 0 - 0; w, = 1/4, i / 0; e, = 1/4, i ± n/2; e n/2 = 0. Then we have (cf. (5.2.5)): 

R(i,-1) = w u R(i, 0)=l/2, R(i, 1) = e, (5.2.7) 

or 

= l w i */ 2 e >] (5.2.8) 

We can also write [i?] = [1 2 l]/4 and stipulate that stencil elements that refer to values 

of u at points outside G are to be replaced by 0. 

The relation between the stencil of an operator and that of its adjoint 

For prolongation operators, a nice definition of stencil notation is less obvious than for restric- 
tion operators. As a preparation for the introduction of a suitable definition we first discuss 

the relation between the stencil of an operator and its adjoint. Define the inner product on 
U in the usual way: 

(u,v)= u i v i (5.2.9) 

where u and v are defined to be zero outside G. Define the transpose A* of A : u -* U i 
the usual way by 

(Au,v) - (u, A*v), Vu,t?ef/ 

Defining A(i,j) = 0 for i £ G or j / S A we can write 


in 


(5.2.10) 


(Au,v) = '£ E A(i,j)u i+jV> = E E A(i,k-i)u k v, 


i t jeZ 




E u k E A(i,k - i)vi = (u,A*v) 


with 


(A v) k - ^2 A(i,k-i) Vi - ^2 A(i + k, -i)v k+ , = ^ A*(k,i)v k+i 


ieZ 


,eZ a 


(5.2.11) 


(5.2.12) 


ie2Z d 


Hence, we obtain the following relation between the stencils of A and A*: 

A*(k,i) = A(k + i, -i) (5.2.13) 

Stencil notation for prolongation operators 

If R U U, then R* : U -» U is a prolongation. The stencil of R * is obtained in similar 
fashion as that of A*. Defining R(i,j) = 0 for i ft G or j £ S R , we have 


(Ru,v) = J2 E R(i,j)u 2t+J Vi = E E R(i,k-2i)u k Vi 
i,j€/Z d i,keZ d 

= E «* E R(i,k-2i)v { = (u,R*v) 

keZ d teZ d 


(5.2.14) 
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with R* :U -*• U defined by 


(5.2.15) 


(. R*v) k = R(^ k ~ 2i ) d i 

i eZ d 

Equation (5.2.15) shows how to define the stencil of a prolongation operator P:U -+ U: 

(Pu)i= -***Ci* * — 2i)“i (5.2.16) 

jeZ d 

Hence, a convenient way to define P is by specifying P*. Equation (5.2.16) is the desired 

stencil notation for prolongation operators. * , ■ „j 

Suppose a rule has been specified to determine Pu for given u , then P (k, m) can be obtain 
as follows. Choose u = 6 k as follows 


(5.2.17) 

(5.2.18) 


4 = 1, 4 = 0, j^k 
Then (5.2.16) gives P*(k, i-2k) = ( P6)i , or 

P*(k,j) = {Pfi k )2k+ji k € G, j G G. 

In other words, [P*]fc is precisely the image of 6 k under P. 

The usefulness of stencil notation will become increasingly clear in what follows. 

Exercise 5.2.1 Verify that (5.2.13) and (5.2.15) imply that, if A and R are represented by 
matrices, A* and iT follow from A and R by interchanging rows and columns. (Remark: for 
d = 1 this is easy; for d > 1 this exercise is a bit technical in the case of R). 

Exercise 5.2.2 Show that if the matrix representation of A : U - U is symmetric, then its 
stencil has the property A(k,i ) = A(k + i, — *)• 

5.3 Interpolating transfer operators 

We begin by giving a number of examples of prolongation operators, based on interpolation. 

Let d = 1, and let G and G be vertex-centred (cf. Figure 5.1.1). Defining P : U -> U by 
linear interpolation, we have 


(Pu)ii = Ui, (Pu) 2 i+1 - + u i+ 1) 


Using (5.3.1) we find that the stencil of P is given by 


[P*) = \[12\] 


(5.3.1) 


(5.3.2) 
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In two dimensions, linear interpolation is exact for functions /(*„x 2 ) = l,x„x„ and takes 
p e in triangles, cf. Figure 5.3.1. Choosing triangles ABD and ACD for interpolation, one 
ains u A - u A , u a - ~(u A + u B ), u e = \(u A + u B ) etc. Alternatively, one may choose 


C d D 

bee 

A a B 


Figure 5.3.1: Interpolation in two dimensions, vertex-centered grids. (Coarse grid point- 
capital letters; fine grid points: capital and lower case letters.) 

triangles ABC and BDC, which makes no essentia] difference. Bilinear interpolation is exact 
or unc ions f(x i,ar 2 ) - and takes place in the rectangle ABCD The only 

difference with linear interpolation is that now = \(u A + u B + u c + u D ). In other words: 

«2t+e,+ e2 - 4 {U, + U t+ei + «, + e 2 -f «,+ e ,+e 2 ), with e, = (1,0) and e 2 = (0, 1) 
lhe stencil for bilinear interpolation is 


i 1 2 1 

1^1 = 4 2 4 2 ( 5 . 3 . 3 ) 

[ 1 2 1 

For the three-dimensional case, see [141]. 

Restrictions 

We can be brief about restrictions. One may simply take 

R = aP * ( 5 . 3 . 4 ) 

with <7 a suitable scaling factor. The scaling of R, i.e. the value of j s important. 

If Ru is to be a coarse grid approximation of u (this situation occurs in non-linear multigrid 
methods, which will be discussed later, then one should obviously have T ■ R(i i) = 1 If 
however, R is used to transfer the residual r to the coarse grid, then the Correct value of 
R(C Jj I depends on the scafing of the coarse and fine grid problems. The rule is that the 
coarse grid problem should be consistent with the differential problem in the same way as 

the fine grid problem. This means the following. Let the differential equation to be solved be 
denoted as 


(5.3.3) 


and the discrete approximation on the fine grid by 


(5.3.5) 


(5.3.6) 
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Suppose that (5.3.6) is scaled such that it is consistent with h a Lu = with ha measure of 
the mesh-size of G. Finite volume discretization leads naturally to o - tn ti d the number 
of dimensions: often (5.3.6) is scaled in order to get nd of div, suras by h. Let the discrete 
approximation of (5.3.5) on the coarse grid G be denoted by 


Au = Rb 


(5.3.7) 


and let A approximate h a L. Then Rb should approximate h a s. Since b approximates h a s, 
we find a scaling rule, as foDows. 


Rule scaling of R: 


: ER{iJ) = (h/h) a 


(5.3.8) 


We emphasize that this rule applies only if R is to be applied to right-hand sides and/or 
residuals. Depending on the way the boundary conditions are implemented ^ boundaries 
a may be different from the interior. Hence the scaling of R should be different at the 
boundary. Another reason why £, R(iJ) may come out different at the boundary is that 
use is made of the fact that due to the boundary conditions the residual to be restricted 

known to be zero in certain points. ...... 

A restriction that cannot be obtained by (5.3.4) with interpolating prolongation is injection. 


( Ru)i = ou 2 i 


(5.3.9) 


Accuracy condition for transfer operators 

The proofs of mesh-size independent rate of convergence of MG assume that P and R sa is y 
certain conditions [21], [57]. The last reference (p. 149) gives the following simple condition. 


mp + nip > 2m 


(5.3.10) 


A necessary condition (not discussed here) is given in [66]. Here orders m P ,mpof P an 
are defined as the highest degree plus one of the polynomials that are interpolated exactly by 
P or siT, respectively, with « a scaling factor that can be chosen freely, and 2m ls the or J 
of the partial differential equation to be solved. For example, (5.3.9) has mp - 0, (5.3.3) has 
m P = 2. Practical experience (see e.g. [139]) confirms that (5.3.10) is necessary. 
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Operator-dependent transfer operators 

If the coefficients m the differentia] equations are discontinuous across certain interfaces be- 
tween subdomains of different physical properties, then u £ C^fl), and linear interpolation 
across discontinuities in u a is inaccurate. (See [141] for more details). Instead of interpola- 
tion, operator-dependent prolongation has to be used. Such prolongations aim to approximate 
the correct jump condition by using information from the discrete operator. They are required 

j- Verte * " centered multigrid, but not in cell-centered multigrid, as shown in [141] where 
a lull discussion of operator-dependent transfer operators may be found. 

5.4 Coarse grid Galerkin approximation 

The problem to be solved on the fine grid is denoted by 

Au = f (5.4.1) 

The two-grid algorithm (2.3.14) requires an approximation A of A on the coarse grid. There 
are basically two ways to chose A, as already discussed in Chapter 2. 

(l) Discretization coarse grid approximation DC A): Uke A, A is obtained by discretization 
of the partial differentia] equation. 

(ii) Galerkin coarse grid approximation ( GCA ): 

A = RAP (5.4.2) 

A discussion of (5.4.2) has been given in Chapter 2. 

The construction of A with DCA does not need to be discussed further. We will use stencil 
notation to obtain simple formulae to compute A with GCA. The two methods will be 
compared, and some theoretical back-ground will be given. 

Explicit formula for coarse grid operator 

The matrices R and P are very sparse and have a rather irregular sparsity pattern. Stencil 
notation provides a very simple and convenient storage scheme. Storage rather than repeated 
evaluation is to be recommended if R and P are operator-dependent. We will derive formulae 
tor A using stencil notation. We have (cf. (5.2.16)) 

(Pu) = p *(j, i ~ 2j)tZj (5.4.3) 

3 

Unless indicated otherwise, summation takes place over Z d . Equation (5.2.1) gives 

( APu )i = ^(*, fc)(P«), + fc = Mh k)P*(j, i + k- 2 j)v.j (5.4.4) 

^ lc 1 
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Finally, equation (5.2.5) gives 

(. RAPu)i = E H(*,m)(APtt) 2 i+m 


= E E E »n)A(2* + TO, k)P*{j,2i + TO + A; - 2j)uj 

m k j 


(5.4.5) 


With the change of variables j = i + n this becomes 

( Au)i = EEE P (*» m) A(2* + TO, k)P*(i +n,m + k- 2 ,n)u i+n (5.4.6) 


m k n 

from which it follows that 


A(i, n) = Y1 £ m ) A ( 2i + m ’ fc ) F *^ + ”» + k - 2n ) ^ 5 - 4 ' 7) 


m k 


For calculation of A by computer the ranges of m and k have to be finite. S A is the structure 
of A as defined in (5.2.2), and S R is the structure R, i.e. 

S R = {j € ZZ d : 3* G G with R{i, j) / 0} (5.4.8) 

Equation (5.4.7) is equivalent to 

A(i,n)= £ H(*, m)A(2i + to, k)P*(i + n, to + fc — 2n) (5.4.9) 

m eS R keS A 

With this formula, computation of A is straightforward, as we will now show. 

Calculation of coarse grid operator by computer 

For efficient computation of A it is useful to first determine S A - This can be done with the 
following algorithm 


Algorithm STRURAP 

comment Calculation of 5^ 
begin 5^ =/0 

for q G Sp * do 

for m G Sp do 

for k G Sj | do 

begin n = ( m + k — <?)/2 

if (n G then S^ = S^Un 

end 


od od od 
end STRURAP 
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Having determined S A it is a simple matter to compute A. This can be done with the fol- 
lowing algorithm. 


Algorithm CALRAP 

comment Calculation of A 
begin A = 0 

/or n 6 S d° 

for m G Sp do 

for k G S £ do 

q = m + k - 2n 
if q G Sp* then 

G\ = {i 6 G : 2i + m e G} n {i 6 G : i + n £ G) 
for i e G\ do 

A(i, n) = A(i, n) + R(i, m)A(2i + m, k)P*(i + n, q) 

od od od 
end CALRAP 

Keeping computation on vector and parallel machines in mind, the algorithm has been de- 
signed such that the innermost loop is the longest. 

To illustrate how (7i is obtained we given an example in two dimensions. Let G and G be 
given by 

G = {i G Z 2 : 0 < i\ < 2 n 1? 0 < 12 ^ 2n 2 } 

G = {* 6 Z 2 : 0 < »i < »,, 0 < *2 < n 2 } 

Then i E G\ is equivalent so 

max(-Ja,-m a / 2,0) < i a < min(n a - m a /2,n a - j a ,n a ) a = 1,2 
It is easy to see that the inner loop vectorizes along grid lines. 

Compansion of discretization and Galerkin coarse grid approximation 

Although DCA seems more straightford, GCA has some advantages. The coarsest grids em- 
ployed in multigrid methods may be very coarse. On such very coarse grids DCA may be 
unreliable if the coefficients are variable, because these coefficients are sampled in very few 
points. An example where multigrid fails because of this effect is given in [137]. The situation 
can be remedied by not sampling the coefficients pointwise on the coarse grids, but taking 
suitable averages. This is, however, precisely that GCA does accurately and automaticaly. 
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For the same reason GCA is to be nsed for interface problems (discontinuous coefficients), 
in which case the danger of pointwise sampling of coefficients is most obvious. Another a - 
vantage of GCA is that it is purely algebraic in nature; no use in made of the under ymg 
differential equation. This opens the possibility of developing autonomous or black ox 
multigrid subroutines, requiring as input only a matrix and right-hand side. On the other 
hand, for non-linear problems and for systems of differential equations there is no general way 
to implement GCA. Both DCA and GCA are in widespread use. 

Structure of coarse grid operator stencil . c 

Galerkin coarse grid approximation will be useful only if S A is not (much) larger than a, 
otherwise the important property of MG, that computing work \sjTovorUond to he number 
of unknowns, may get lost. For examples and further discussion of CGA, including the possible 
loss of the K- matrix property on coarse grids, see [141]. 


6 Multigrid algorithms 

6.1 Introduction 

The order in which the grids are visited is called the multigrid schedule Several schedules will 
be discussed. All multigrid algorithms are variants of what may be called the basic mutignd 
algorithm. This basic algorithm is nonlinear, and contains linear multignd as a special case. 
The most elegant description of the basic multigrid algorithm is by means of a recursive 
formulation. FORTRAN does not allow recursion, thus we also present a non-recursive ior- 
mulation. This can be done in many ways, and various flow diagrams have been presented m 
the literature. If, however, one constructs a structure diagram not many possibilities remain, 
and a well structured non-recursive algorithm containing only one goto statement resu s. 
The decision whether to go a finer or to a coarser grid is taken in one place only. 

6.2 The basic two-grid algorithm 

Preliminaries . k , , 

Let a sequence {G k : k = 1,2,..., A} of increasingly finer grids be given. Let U be the 

set of grid functions G k — K on G k \ a grid function U £ U stands for m functions in 
the case where we want to solve a set of equations for m unknowns Let there be given 
transfer operators P k : U k ~ l -*• U k (prolongation) and R :U -* U (restriction). Let 
the problem to be solved on G k be denoted by 


L k (u k ) = b k 


( 6 . 2 . 1 ) 
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The operator L may be linear or non-linear. Let on every grid a smoothing algorithm be 
defined, denoted by S{u,v, f ,v,k). S changes an initial guess u k into an improved approxi- 
mation v with right-hand side f k by v k iterations with a suitable smoothing method. The 
use of the same symbol u k for the solution of (6.2.1) and for approximations of this solution 
will not cause confusion; the meaning of u k will be clear from the context. On the coarse 
grid G 1 we sometimes wish to solve (6.2.1) exactly; in general we do not wish to be specific 
about this, and we write S(u,v, /, •, 1) for smoothing or solving on G 1 . 

The nonlinear two-grid algorithm 

Let us first assume that we have only two grids G k and G k ~\ The following algorithm is a 
generalization of the linear two-grid algorithm discussed in Section 2.3. Let some approxima- 
tion u of the solution on G k be given. How u k may be obtained will be discussed later. The 
non-linear two-grid algorithm is defined as follows. Let f k = b k . 

Subroutine TG (u^u^f, k) 
comment nonlinear two-grid algorithm 
begin 

S(u,u,f y v, k) 
r k - f k - L k (u k ) 

Choose 

f k - 1 = + 

S(u,u,f,*,k- 1) 

u k = u k + (l/s k ^)P k (u k - 1 - u k -') 

S(u, u , f,n,k) 
end of TG 

A call of TG gives us one two-grid iteration. The following program performs ntg two-grid 
iterations: 

Choose u k 
f k = b k 

for i = 1 step 1 until ntg do 

TG ( u,u,f,k ) 

u = u 

od 

Discussion 

Subroutine TG is a straightforward implementation of the basic multigrid principles discussed 
in Chapter 2, but there are a few subtleties involved. 

We proceed with a discussion of subroutine TG. Statement (1) represents i/ k smoothing it- 


( 1 ) 

( 2 ) 

( 3 ) 

( 4 ) 

( 5 ) 

( 6 ) 
( 7 ) 
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erations (pre-smoothing), starting from an initial gness u k . In (2) the residual r* is com- 
puted' r k is going to steer the coarse grid correction. Because short waveleng accu y 
already achieved fn u k mas. not get lost, u k is to be kept, and a correct, on (containing 
■long wavelength information') is to be added to u k In the non-linear case, r cannot be 
taken for the right-hand side of the problem for 6u k ; L({u k ) = r‘ might not even have a 
solution. For the same reason, R k ~' r k cannot be the right-hand side for the coarse gr 
problem on G k ~' . Instead, it is added in (4) to with u " an approsrmatmn 

to the solution of (6.2.1) in some sense (e.g. Pu k ~' ~ solution of equation (6.2.1). Ob- 
viously, L k -'(u 1 -') = has a solution, and if R r is not too large, en 

_ L k ~ l (i i k ~ l ) + R k ~ x r k can also be solved, which is done in statement (5) (ex- 
actly or approximately). 

R k ~'r k will be small when u k is close to the solution of equation (6.2.1), i.e. when the 
algorithm is close to convergence. In order to cope with situations where R r is not smal 
enough, the parameter s k . x is introduced. By choosing s k x small enough one can bring / 

arbitrarily close to L k ~ l {u k ~ 1 )- Hence, solvability of L "(«') = / ca « * ^k-l t 
Furthermore, in bifurcation problems, u k ~' can be kept on the same branch u by 
means of s k j . In (6) the coarse grid correction is added to u . Omission of the factor / fc-i 
luld meaiT that only part of the coarse grid correction is added to 

damping of the coarse grid correction; this would slow down convergence. Finally, statement 
(7) represents fi k smoothing iterations (post-smoothmg). 

The linear two-grid algorithm , , . , 

It is instructive to see what happens when L k is linear. It is reasona e ^ o assume a ’ e 

L k ~ x is also linear. Furthermore, let us assume that the smoothing method is linear, that is 

to say, statement (5) is equivalent to 


u 


k - 1 


= u 


k - 1 


fc-1 / fk - 1 


+ B k -\f 


— L k i u k 1 


) 


with some linear operator. With f k 1 from statement (4) this gives 


u k 1 = u 


+ s k -\B k ~ x R k ~ l T k 


( 6 . 2 . 2 ) 


(6.2.3) 


Statement (6) gives 


u' 


= u k + P k B k ~ l R‘ 


i k d A: — 1 1 


(6.2.4) 


and we see that the coarse grid correction P k B k " R k -'r k is independent ofttechoice of 
s k -j and u k ~ 1 in the linear cas. Hence, we may as well choose s k . a = 1 and u - U in tne 
linear case. This gives us the following linear two-grid algorithm. 
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Subroutine LTG (u y u,f y k) 
comment linear two-grid algorithm 
begin 

S(u, u,f,v,k) 

r k = f k - L k u k 
fk-i _ R k-\ r k 

•u fc_1 = 0 

S(u,u, 1) 

u k = u k + p k u k - 1 

end of LTG 

Choice of u k ~ 1 and s k -\ 

There are several possibilities for the choice of -u* -1 . One possibility is 

u*- 1 = R k ~\ k (6.2.5) 

where R is a restriction operator which may or may not be the same as R k ~ l . 

With the choice = 1 this gives us the first non-linear multigrid algorithm that has 
appeared, the FAS (full approximation storage) algorithm proposed by Brandt [20], The more 
general algorithm embodied in subroutine TG, containing the parameter s k -\ and leaving the 
choice of open, has been proposed by Hackbusch [54], [49], [57]. In principle it is possible 
to keep u k -i fixed, provided it is sufficiently close to the solution of L k ~ l (u k ~ x ) = &* _1 . This 
decreases the cost per iteration, since L k 1 (-u* needs to be evaluated only once, but the 
rate of convergence may be slower than with tt* -1 defined by (5). We will not discuss this 
variant. Another choice of u A_1 is provided by nested iteration, which will be discussed later. 
Hackbusch [54], [49], [57] gives the following guidelines for the choice of ii* -1 and the param- 
eter s k i. Let^the non-linear equation L k ~ l (u k - 1 ) = f k ~ x be solvable for ||/ A:-1 || < p k 
Let ||Z 1 ( wA_1 )|| < Pk- 1/2. Choose such that ||« fc -i || < p k - 1/2, for example: 

=\pk-xl\\R k - X r k \\. ( 6 .2.6) 

Then \\f k -'\\ < pk- 1 , so that the coarse grid problems has a solution. 

6.3 The basic multigrid algorithm 
The recursive non-linear multigrid algorithm 

The basic multigrid algorithm follows from the two-grid algorithm by replacing the coarse 
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grid solution statement (statement ( 5 ) in subroutine TG) by 7* multigrid iterations. This 
leads to 

Subroutine MG 1 (u,u, 

comment recursive non-linear multigrid algorithm 

begin 

if (k eq 1) then 
S(u,u,f,»,k) 

else 

S(u,u, f , V) fc) 
r k = f k - L k (u k ) 

Choose ii k ~ 1 ,Sk-i 
f k - 1 = L k ~'{u k ~ x ) + s k -\R k - l r k 
for i — 1 step 1 until 7* do 
MG 1 (u,u,f,k- 1,7) 
od 

u k = u k + (l/s fc _ 1 )P* : (u ' c - 1 - 'U fc ' 1 ) 

S (ti, n, /, /x, fc) 
endif 
end of MG1 

After our discussion of the two-grid algorithm, this algorithm is self explanatory. ^ 

The following program carries out nmg multigrid iterations, starting on the finest grid G . 

Program 1 : 

Choose u K 

f K = b K 

for i = 1 step 1 until nmg do 
MG 1 (u,u,f,K, 7) 
od 


( 1 ) 

( 2 ) 

( 3 ) 

( 4 ) 

( 5 ) 

( 6 ) 

( 7 ) 

( 8 ) 
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The recursive linear multigrid algorithm 

The linear multigrid algorithm follows easily from the linear two-grid algorithm LTG: 

Subroutine LMG ( u,u,f,k ) 

comment recursive linear multigrid algorithm 
begin 

if ( k — 1) then 
S(u,u, /, »,k) 

else 

S(u,u, f,v,k) 

r k = f k - L k u k 
fk-\ _ R k-\ r k 

= 0 

fori- 1 step 1 until 7 k do 
LMG (u,u,f,k - 1) 

u k ~ 1 = U k ~ 1 
od 

u k = u k + P k u k ~ 1 

S(u,u,f,p,k) 

endif 
end LMG 


Here u plays the role of an initial guess. 

Multigrid schedules 

The order in which the grids are visited is called the multigrid schedule or multigrid cycle. 
f the parameters 7 k, k- 1 , 2 ,..., K - 1 are fixed in advance we have a fixed schedule ■ if 
depends on intermediate computational results we have an adaptive schedule. Figure 6.3.1 
shows the order in which the grids are visited with 7fc = 1 and 7* = 2, k = 1, 2, ..., K- 1 , in the 
case K = 4. A dot represents a smoothing operation. Because of the shape of these diagrams 
these schedules are called the V-, W- and sawtooth cycles, respectively. The sawtooth cycle is 
a special case of the V-cycle, in which smoothing before coarse grid correction (pre-smoothing) 
is deleted. A schedule intermediate between these two cycles is the F-cycle. In this cycle coarse 
grid correction takes place by means of one F-cycle followed by one V-cycle. Figure 6 3 2 gives 
a diagram for the F-cycle, with K - 5. & 


Recursive algorithm for V-, F- and W-cycle 

A version of subroutine MG1 for the V-, W- and F-cycles is as follows. The parameter 7 is 
now an integer instead of an integer array. 


k 



Subroutine MG2 (u,u, f 7) 

comment nonlinear multigrid algorithm V-, W-, or F-cyc e 
begin 

if (k eq 1 ) then 

S(u,u,f,*,k) 

if (cycle eq F) then 7 = 1 endtf 

else 

A 

for i = 1 step 1 until 7 do 

MG2 1,7) 

od 

g 

if (keq K and cycle eq F ) then 7 = 2 endif 

endif 

end MG2 

Here rl aad B represent statements (2) to (5) and (7) and (8) in snbrontine MG1. 
following program carries out nmg V-, W-, or F-cycles. 


The 
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Figure 6.3.2: F-cycle diagram. 


Program 2 : 

Choose u h 
f K = b K 

if (cycle eq W or cycle eq F) then 7 = 2 else 7 = 1 
for i = 1 step 1 until rung do 

MG2 (u,u,/, A , 7) 

od 

Adaptive schedule 

An example of an adaptive strategy is the following. Suppose we do not carry out a fixed 
number of multigrid iterations on level G^ y but wish to continue to carry out multigrid 
interactions, until the problem on G ^ is solved to within a specified accuracy. Let the accuracy 
requirement be 

\\L k (u k ) - /*|| <£* = 6 k \\L k+ '(u k+l ) - f k +'\\ ( 6 . 3 . 1 ) 

with 6 6 ( 0 , 1 ) a parameter. 

At first sight, a more natural definition of e k would seem to be e k = £||/ fe ||. Since f k does 
not, however, go to zero on convergence, this would lead to skipping of coarse grid correction 
when u k+1 approaches convergence. Analysis of the linear case leads naturally to condition 
(6.3.1). An adaptive multigrid schedule with criterion (6.3.1) is implemented in the following 
algorithm. In order to make the algorithm finite, the maximum number of multigrid iterations 
allowed is 7 . 
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Subroutine MG3 /,fc) 

comment recursive nonlinear multigrid algorithm with adaptive 
schedule 

begin 

if ( k eq 1) then 

S (tlr, U, /, ^ 

else 

A 

(1) 4-i = ll* , *ll- £ * 

e k -\ = 
rik - 1 = 7 

while (4_i > 0 and nk-i > 0) 

MG3 (u,u,f,k- 1) 

7lfc_i = 71fc_i — 1 

4-i = \\L k ~ 1 (u k ~ 1 )-f k ~ 1 \\ - £k-i 

od 

B 

endif 

end MG3 

Here A and B stand for the same groups of statements as in subroutine MG2. The purpose 
of statement (1) is to allow the possibility that the required accuracy is already reached by 
pre-smoothing on G k , so that coarse grid correction can be skipped. The following program 
solves the problem on G K within a specified tolerance, using the adaptive subroutine MG3: 

Program 3 : 

Choose u K 

f K = b K ; = tol * ||6 k || ; t K = \\L K (u K ) - b K \\ - e k 

n — nmg 

while (tx > 0 and n > 0) do 
MG3 (u,uJ,K) 

r tK = \\L K (u K )-b K \\-£ K 

od 

The number of iterations is limited by mng. 

Storage requirements 

Let the finest grid G K be either of the vertex-centered type given by (5.1.1) or of the cell- 
centered type given by (5.1.2). Let in both cases n a = n« ^ = m a • 2 . Let the coarse 
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grids G k k - K - 1, K - 2, 1 be constructed by successive doubling of the mesh-sizes h a 
( standard coarsening). Hence, the number of grid-points N k of G k is 


d 

N k= + M2 kd 

01 = 1 


(6.3.2) 


in the vertex-centered case, with 


and 


d 

M = n rn a , 

Ct = 1 


N k = M2 kd 

in the cell-centered case. In order to be able to solve efficiency 
desirable that m a is small. Henceforth, we will not distinguish 
and cell-centered case, and assume that N k is given by (6.3.3). 


(6.3.3) 

on the coarsest grid G 1 it is 
between the vertex-centered 


It is to be expected that the amount of storage reqired for the computations that take 

place on G is given by Cl N k , with c, some constant independent of k. Then the total amount 
of storage required is given by 


K 2 d 
ci 12 N * - » , c i N k 

k=\ z 1 


(6.3.4) 


Hence, as compared to single grid solutions on method selected, the use of multigrid increases 
the storage required by a factor of 2 d /(2 d -l), which is 4/3 in two and 8/7 in three dimensions, 
so that the additional storage requirement posed by multigrid seems modest. 

Next, suppose that semi-coarsening (cf. Section 7.3) is used for the construction of the coarse 

grids G ,k < K. Assume that in one coordinate direction the mesh-size is the same on all 
grids. Then 


N k = M 2 a '+^-i) 

and the total amount of storage required is given by 


(6.3.5) 


K 


Cl £ N k ~ 

Jt=l 


2^-1 

2^-1 _ j Cl Nk 


(6.3.6) 


Now the total amount of storage required by multigrid compared with single grid solution on 
increases by a factor 2 in two and 4/3 in three dimensions. Hence, in two dimensions the 
storage cost associated with semi-coarsening multigrid is not negligible. 


99 


Computational work , 

We wffl estimate the computational work of one iteration with the fixed schedule algori hm 
MG2 A close approximation of the computational work w k to be performed on G k will be 
w , = C2 N k , assuming the number of pre- and post-smoothings v k and /x fc are independent ol 
k and that the operators L k are of similar complexity (for example, in the linear case, L are 
matrices of equal sparsity). More precisely, let us define w k to be all computing work involved 
in MG2 (it, u, /, k ), except the recursive call of MG2. Let be all work involved in MG2 
( u,u,f,k ). Let 7fc = 7, fc = 2,3,..., A - 1, in subroutine MG2 (e.g., the V- or W-cycles). 

Assume standard coarsering. Then 



W k = c 2 M 2 kd + ^W k -\ 

(6.3.7) 

In [141] it is shown that if 

7 = 7/2 d < I 

(6.3.8) 

then 

W K < W = 1/(1 -7) 

(6.3.9) 


with Wk = W k /(c 2 N k ) ■ , 

The following conclusions may be drawn from (6.3.10). W K is the ration of multigrid wor 
work on the finest grid. The bulk of the work on the finest grid usually consists of smoothing. 
Hence, W K - 1 is a measure of the additional work required to accelerate smoothing on the 
finest grid G K by means of multigrid. 

If 7 > 1 the work Wk is superlinear in the number of unknowns N K , see [141]. 


If 7 < 1 equation (6.3.9) gives 

W k <c 2 N k /( 1-7) (6.3.10) 

so that W K is linear in N K - It is furthermore significant that the constant of proportionality 
C2 /(i _ 7) is small. This because c 2 is just a little greater than the work per grid point of the 
smoothing method, which is supposed to be a simple iterative method (if not, multigrid 

is not applied in an appropriate way). Since an (perhaps the main) attractive feaure of 
multigrid is the possibility to realize linear computational complexity with small constant ot 
proportionality, one chooses 7 < 1, or 7 < 2 d . In practice it is usually found that 7 > 2 does 
not result in significantly faster convergence. The rapid growth of W K with 7 means that it 
is advantageous to choose 7 < 2, which is why the V- and W-cycles are widely used. 

The computational cost of the F-cycle may be estimated as follows. In Figure 6.3.3 the 
diagram of the F-cycle has been redrawn, distinguishing between the work that is done on 
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G preceding coarse grid correction (pre-work, statements A in subroutine MG2) and after 
coarse grid correction (post-work, statements B in subroutine MG2). The amount of pre- 
and post-work together is c 2 M2 kd , as before. It follows from the diagram, that on G k the 
cost of pre- and post-work is incurred jk times, with jk = K - k + I, k = 2, 3,..., A', and 

to = K .~ h F ° r convenience we redefine j x = K, bearing our earlier remarks on the inaccuracy 
and unimportance of the estimate of the work in G 1 in mind. One obtains 


K 


Wk = c 2 M YXK — k + \)2 kd 
k= 1 


We have 


A ^ 2l K+ W J od 

as is checked easily. It follows that 

W K = c 2 M( 2 d ( K+ V + K + 1 - K2 d )j(2 d - if 


(6.3.11) 


(6.3.12) 


k 

5 

4 

3 

2 

1 



Figure 6.3.3: F-cycle (o pre-work, • post-work). 

so that 

W K < W = 1/(1 - 2~ d f (6.3.13) 

Table 6.3.1 gives W as given by (6.3.9) and (6.3.13) for a number of cases. The ratio of 
multigrid over single grid work is seen to be not large, especially in three dimensions. The 
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d 

2 

3 

V- cycle 

4/3 

8/7 

F-cycle 

16/9 

64/49 

W-cycle 

2 

4/3 

7 = 3 

4 

8/5 


Table 6.3.1: Values of W, standard coarsening 

F-cycle is not much cheaper than the W-cycle. In three dimensions the cost of the V-, F- and 
W-cycles is almost the same. 

Suppose next that semi-coarsening is used. Assume that in one coordinate direction the 
mesh-size is the same on all grids. The number of grid-points N k of G is given by (6.3.5). 
With 7 *. = 7 , k = 2 , 3, ..., K - 1 we obtain 

W k = c 2 M2 K+k(d ~ 1) +jW k -\ (6.3.14) 

Hence W k is given by (6.3.8) and W by (6.3.9) with 7 = l/2 d ~ l . For the F-cycle we obtain 

W K = c 2 M2* - k + l) 2 fc ( d - 1 > (6.3.15) 

k—\ 


Hence 


W K < W = 1/(1 -2 l ~ d ) 2 


d 

2 

3 

V- cycle 

2 

4/3 

F-cycle 

4 

16/9 

W-cycle 

- 

2 

7=3 

- 

4 


Table 6.3.2: Values of W , semi-coarsening 

Table 6.3.2 gives W for a number of cases. In two dimensions 7 = 2 or 3 is not useful, because 
7 > 1. It may happen that the rate of convergence of the V-cycle is not independent of the 
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mesh-size, for example if a singular perturbation problem is being solved (e.g. convection- 
lffusion problem with £ < 1), or when the solution contains singularities. With the W-cycle 
we have 7 = 1 with semi-coarsening, hence W k = K. In practice, K is usuaUy not greater 
than 6 or 7, so that the W-cycle is still affordable. The F-cycle may be more efficient 


Work units 

The ideal computing method to approximate the behaviour of a given physical problem in- 
volves an amount of computing work that is proportional to the number and size of the 
physical changes that are modeled. This has been put forward as the ’golden rule ofcompu- 

AS haS been em P hasized ^ Brand ^ a "umber of publications, e.g. 
l^UJ, [21J, [22J, [16J, this involves not only the choice of methods to solve (6.2.1), but also the 
choice of the mathematical model and its discretization. The discretization and solution pro- 
cesses should be inter wined, leading to adaptive digitization. We shall not discuss adaptive 
methods here, but regard (6.2.1) as given. A practical measure of the minimum computing 
work to solve (6.2.1) is as follows. Let us define one work unit (WU) as the amount of comput- 
rng work required to evaluate the residual L K (u K ) - b K of Equation 6.2.1) on the finest grid 
. hen it is to be expected that (6.2.1) cannot be solved at a cost less than few WU, and 
one should be content if this is realized. Many publications show that this goal can indeed 
be achieved with multigrid for significant physical problems, for example in computational 
ui dynamics. In practice the work involved in smoothing is by far the dominant part of the 
total work. One may, therefore, also define one work unit, following [20], as the work involved 

in one smoothing iteratio " the finest grid G K . This agrees more or less with the first 

definition only if the smoothing algorithm is simple and cheap. As was already mentioned, if 

t 1S J* .T* th ! , CaS !^ U i! igrid ' S n0t ap P lied in an appropriate way. One smoothing iteration 
on G then adds 2 » WU to the total work. It is a good habit, followed by many authors, 

to pubbsh convergence histones in terms of work units. This facilitaties comparisons between 
methods, and helps in developing and improving multigrid codes. 


6.4 Nested iteration 
The algorithm 

Nested iteration , also called full multigrid (FMG, [22], [16]) is based on the following idea 
When no a prion information about the solution is available to assist in the choice of the 
initial guess u on the finest grid G K , it is obviously wasteful to start the computation on 
the finest grid, as is done by subroutines MGi, i = 1,2,3 of the preceding section. With 
an unfortunate choice of the initial u h , the algorithm might even diverge for a nonlinear 
problem. Computing on the coarse grids is so much cheaper, thus it is better to use the 
coarse grids to provide an informed guess for u h . At the same time, this gives us a choice 
or u ’ < A . Nested iteration is defined by the following algorithm. 
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( 1 ) 

( 2 ) 

( 3 ) 


Program 1 

comment nested iteration algorithm 
Choose u 1 

S (n, n, f i 1 

for k = 2 step 1 until K do 

k ~ k 1 

U = U = jT U 

for i = 1 step 1 until 7* do 
MG (u,u,f,k) 

od 

od 


Of course, the value of 7fc inside MG may be different from 7 fc- 


Choice of prolongation operator 

The prolongation operator P* does not need to be identical to P . In fact there may be 
good reason to choose it differently. As discussed in [57] (for a simplified analysis see [141]), 

it is often advisable to choose P such that 

mp > m c (6.4.1) 

where mp is the order of the prolongation operator as defined in Section 5.3, and m is the 
order of consistency of the discretizations L k , here assumed to be^the same on all grids. Of- 
ten m c = 2 (second-order schemes). Then (6.4.1) implies that P is exact for second-order 
polynomials. 

Note that nested iteration provides u k \ this is an alternative to (6.2.5). 

As discussed in [57] and [141], if MG converges well then the nested iteration algorithm 
results in a u K which differs from the solution of (6.2.1) by an amount of the order of the 
truncation error. If one desires, the accuracy of u K may be improved further by following 
the nested iteration algorithm with a few more multigrid iterations. 

Computational cost of nested iteration . , . 

Let = 7 k = 2 3, ..., K, in the nested iteration algorithm, let Wk be the work involved m 

MG(u,u,f,k), and assume for simplicity that the (negligible) work on G l equals W x . Then 
the computational work W ni of the nested iteration algorithm, neglecting the cost of P , is 
given by K 

W ni = £ i w k 


(6.4.2) 
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Assume inside Mg lk = 7, k = 2,3,... ,K and let 7 = 1 /2 d < 1 . Note that 7 and 7 may be 
different. Then it follows from (6.3.9) that 


w ■ 


< 7‘ 




K 


c 2 7 




(1 - t )(1 ~ 2 -^) 




(6.4.3) 


Defining a work unit as 1 WU = c 2 N K , i.e. approximately the work of (v + p) smoothing 
iterations on the finest grid, the cost of a nested iteration is 


^.=7/[(l-7)(l-2i]WU (6.4.4) 

Table 6.4.1 gives the number of work units required for nested iteration for a number of cases. 
The cost of nested iteration is seen to be just a few work units. Hence the fundamental 
property, which makes multigrid methods so attractive: multigrid methods can solve many 
problems to within truncation error at a cost of cN arithmetic operations. Here N is the 
number of unknowns, and c is a constant which depends on the problem and on the multigrid 
method (choice of smoothing method and of the parameters u k ,p k ,^ k ). If the cost of the 
residual b - L (u ) is bN, then c need not be larger than a small multiple of b. Other 
numerical methods for elliptic equations require 0(N a ) operations with a > 1, achieving 
°(N In N) only in special cases (e.g. separable equations). A class of methods which is com- 
petitive with multigrid for linear problems in practice are preconditioned conjugate gradient 
methods. Practice and theory (for special cases) indicate that these require 0(N a ) opera- 
tions, with a = 5/4 in two and a = 9/8 in three dimensions. Comparisons will be given later. 


d 

7 2 

3 

1 16/9 

2 8/3 

64/49 

48/21 


Table 6.4.1: Computational cost of nested iteration in work units; 7 = 1 


6.5 Non-recursive formulation of the basic multigrid algorithm 
Structure diagram for fixed multigrid schedule 

In FORTRAN, resursion is not allowed: a subroutine cannot call itself. The subroutines MG1, 
2, 3 of Section 6.3 cannot, therefore, be implemented directly in FORTRAN. A non-recursive’ 
version will, therefore, be presented. At the same time, we will allow grater flexibility in the 
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decision whether to go to a finer or to a coarses grid. 


Various flow diagrams describing non-recursive multigrid algorithms have been pubbshe , 
for example in [20] and [57]. In order to arrive at a well structured program, we begin y 
presenting a structure diagram. A structure diagram allows much less freedom in the design 
of the control structure of an algorithm than a flow diagram. We found basically only one 
way to represent the multigrid algorithm in a structure diagram ([134], [140]). This structure 
diagram might, therefore, be called the canonical form of the basic multigrid algorithm. The 
structure diagram is given in Figure 6.5.1. This diagram is equivalent to 1 Progr; am 2 calling 
MG2 to do nmg multigrid iterations with finest grid G in Section 6.3. The schedule is fixe 
and includes the V-, W- and F-cycles. Parts A and B are specified after subroutine MG2 in 
Section 6.3. Care has been taken that the program also works as a single grid method for 

K = 1. 


FORTRAN implementation of while clause 

Apart from the while clause, the structure diagram of Figure 6.5.1 can be expressed direct y 
in FORTRAN. A FORTRAN implementation of a while clause is as follows. Suppose we 

have the following program 


while ( n(K ) > 0) do 
Statement 1 
n(K) = ... 

Statement 2 

od 

A FORTRAN version of this program is 


10 if ( n(K ) > 0) then 
Statement 1 
n(K) = ... 
Statement 2 

goto 10 
endif 
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Choose u * and 7 

comment 7 = 1: V-cycle : 7 = 2 ; W-cycIe 
= b * 1 ; k = K ; nx = nmg 
if (cycle eq F ) then 7 = 2 endif 
while (nx > 0) do 

1 


eq 0 or k eq 1 

T 


A 


k eq 1 



k = k — 1 
rik = 7 


S(u,u, f,-,k ) 
if (cycle eq F) then 
7 = 1 

endif 



F 




k = k + 1 

B 

if (cycle eq F) then 
7 = 2 
endif 



3 

II 

9 t- 

1 


Figure 6 . 5 . 1 : Structure diagram of non-recursive multigrid algorithm with fixed schedule 
including V-, W- and F-cycles. 
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The goto statement required for the FORTRAN version of the while clause is the only goto 
needed in the FORTRAN implementation of the structure diagram of Figure 6.5.1. 
FORTRAN implementation is quite obvious, and will not be given. 


Testing of multigrid software . 

A simple way to test whether a multigrid algorithm is functioning properly is to measure the 

residual before and after each smoothing operation, and before and after each visit coarser 
grids If a significant reduction of the size of the residual is not found, then the relevan 
part of the algorithm (smoothing or coarse grid correction) is not functioning P™P orly -f°" 
simple test problems predictions by Fourier smoothing analysis and the contraction number 
of the multigrid method should be correlated. If the coarse grid problem is solved exactly (a 
situation usually approximately realized with the W-cycle) the multigrid contraction number 
should usually be approximately equal to the smoothing factor. 


It°may S however, happen, happen that for a well designed multigrid algorithm the contraction 
number is significant iorse than predicted by the smoothing factor. This may be caused 
by the fact that Fourier smoothing analysis is locally not apphcable. The cause may be a 
local singularity in the solution. This occurs for example when the physical domain has a 
reentrant corner The coordinate mapping from the physical domain onto the computationa 
rectangle is singular at that point. It may well be that the the smoothing method does no 
reduce the residual sufficiently in the neighbourhood of this singularity, a fact that does no 
remain undetected if the testing procedures recommended above are applied The reme > >s 
to apply additional local smoothing in a small number of points in the neighbourhood of the 
singularity. This procedure is recommended in [16], [17], [18], [9], and justifie eore y 
in [110] and [24], This local smoothing is applied only to a small number of points, thus the 
computing work involved is negligible. 


6.6 Remarks on software 

Multigrid software development can be approached in various ways, two of which will be 
examined here. 


The first approach is to develop general building blocks and diagnostic tools, which helps 
users to develop their own software for particular applications without having to start from 
scratch, users will, therefore, need a basic knowledge of multigrid methods. Such software 

tool are described in [26]. 


The second approach is to develop autonomous (black box) programs, for which the user 
has to specify only the problem on the finest grid. A program or subroutine may be called 
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autonomous if it does not require any additional input from the user apart from problem spec- 
ification, consisting of the linear discrete system of equations to be solved and the right-hand 
side. The user does not need to know anything about multigrid methods. The subroutine is 
perceived by the user as if it were just another linear algebra solution method. This approach 
is adopted by the MGD codes [133], [67], [69], [68], [107], [108], which are available in the 
NAG hbrary, and by the MGCS code [36]. 

Of course, it is possible to steer a middle course between the two approaches just outlined 
allowing or requiring the user to specify details about the multigrid method to be used, such 
as offering a selection of smoothing methods, for example. Programs developed in this vein 
are BOXMG [37], [38], [39], the MG00 series of codes [45], [44], [113] which is available 
in ELLPACK [100], MUD PACK [3], [2], and the PLTMG code [11], [10], [12], Exept for 
PLTMG and MGD, the user specifies the linear differential equation to be solved and the 
program generates a finite difference discretization. PLTMG generates adaptive finite ele- 
ment disretizations of non-linear equations, and therefore has a much wider scope then the 
other packages. As a consequence, it is not (meant to be) a solver as fast as the other methods. 


By sacrificing generality for efficiency very fast multigrid methods can be obtained for 
special problems, such as the Poisson or the Helmholtz equation. In MG00 this can be done 
by setting certain parameters. A very fast multigrid code for the Poisson equation has been 
developed in [13], This is probably the fastest two-dimensional Poisson solver in existence. 

If one wants to emulate a linear algebraic systems solver, with only the fine grid matrix and 
right-hand side suplied by the user, then the se of coarse grid Galerkin appoximation (Chap- 
ter 5) is mandatory. Coarse grid Galerkin approximation is also required if the coefficients 
in the differential equations are discontinuous. Coarse grid Galerkin approximation is used 
m MGD, MGCS and BOXMG; the last two codes use operator-dependent transfer operators 
and are applicable to problems with discontinuous coefficients. 


In an autonomous subroutine the method cannot be adapted to the problem, so that user 
expertise is not required. The method must, therefore, be very robust. If one of the smoothers 
that were fund to be robust in Chapter 4 is used, the required degree of robustness is indeed 
obtained for linear problems. 

Non-linear problems may be solved with multigrid codes for linear problems in various 
ways. The problem may be linearized and solved iteratively, for example by Newton method. 
This works well as long as the Jacobian of the non-linear discrete problem is non-singular. 
It may well happen, however, that the given continuous problem has no Frechet derivative. 
In this case the condition of the Jacobian deteriorates as the grid is refined, and the Newton 
method does not converge rapidly or not at all. An example of this situation is given [96], 
[95]. The non-linear multigrid method can be used safely and efficiently, because the global 
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system is not linearized. A systematic way of applying numerical software outside the class 
of problems to which the software is directly applicable is the defect correction approach. In 
[5] and [15] it is pointed out how this ties in with multigrid methods. 

Much software is made available on MGNet. 

6.7 Comparison with conjugate gradient methods 

Although the scope and applicability of multigrid principles are much broader multigrid 
methods can be regarded as very efficient ways to solve linear systems arising 
tion of partial differential equations. As such multigrid can be viewed as a tedmqne to 
accelerate the convergence of basic iterative methods (called smoothers in the multigrid con- 
text) Another powerful technique to accelerate basic iterative methods for linear problems 
that also has come to fruition relatively recently is provided by conjugate gradient an re- 
lated methods. For an introduction to conjugate gradient acceleration of iterative methods, 
see [63], [48] or (for a very brief synopsis) [141]. 

It is surprising that, although the algorithm is much simpler, the rate of convergence of 
conjugate gradient methods is harder to estimate theoretically than for 'mu tignd methods. In 
two^ dimensions, 0{N b ' A ) computational complexity, and probably 0{N ) m three dimen- 

sions seems to hold approximately quite generally for conjugate gradient methods precon >- 
tioned by approximate factorization, which comes close to the 0(N ) of multigrid metho . 

Conjugate gradient acceleration of multigrid . 

The Conjugate gradient method can be used to accelerate any iterative method including 
multigrid methods. If the multigrid algorithm is well designed and fits the problem it wiU 
converge fast, making conjugate gradient acceleration superfluous or even wasteful. If m 
grid does not converge fast one may try to remedy this by improving the algorithm (for 
example introducing additional local smoothing near singularities, or adapting the smoother 
to the problem), but if this is impossible because an autonomous (black box) multigrid code is 
used, or difficult because one cannot identify the cause of the trouble, then conjugate gradient 
acceleration is an easy and often very efficient way out. 

The non-symmetric case . 

A severe limitation of conjugate gradient methods is their restriction to hnear systems with 

symmetric positive definite matrices. A number of conjugate gradient type ™ ethods 
been proposed that are applicable to the non-symmetric case. Although theoretical es- 
timates are available, their rate of convergence is often satisfactory in practice 
methods are CGS (conjugate gradient squared), described in [107], [108], [105] an [ ], an 

GMRES, described in [102], [121], [120], [131], [132]. Good convergence is expected if the 
eigen eigenvalues of A have positive real part, cf. the remarks on convergence in [105] . 
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Comparison of conjugate gradient and multigrid methods 

Realistic estimates of the performance in practice of conjugate gradient and multigrid methods 
by purely theoretical means are possible only for very simple problems. Therefore numerical 
experiments are necessary to obtain insight and confidence in the efficiency and robustness 
of a particular method. Numerical experiments can be used only to rule out methods that 
fad, not to guarantee good performance of a method for problems that have not yet been 
attempted. Nevertheless, one strives to build up confidence by carefully choosing tests prob- 
lems, trying to make them representative for large classes of problems, taking into account 
the nature of the mathematical models that occur in the field of application that one has in 
mind. For the development of conjugate gradient and multigrid methods, in particular the 

subject areas of computational fluid dynamics, petroleum reservoir engineering and neutron 
diffusion are pace-setting. 

Important constant coefficient test problems are (4.5.3) and (4.5.4). Problems with con- 
stant coefficients are thought to be representative of problems with smoothly varying coeffi- 
cients. Of course, in the code to be tested the fact that the coefficients are constant should 
not be exploited. As pointed out in [35], one should keep in mind that for constant coefficient 
problems the spectrum of the matrix resulting from discretization can have very special prop- 
erties, that are not present when the coefficients are variable. Therefore one should also carry 
out tests with variable coefficients, especially with conjugate gradient methods, for which the 
properties of the spectrum are very important. For multigrid methods, constant coefficient 
test problems are often more demanding than variable coefficient problems, because it may 
happen that the smoothing process is not effective for certain combinations of s and (3. This 
fact goes easily unnoticed with variable coefficients, where the unfavourable values of e and 
(3 perhaps occur only in a small part of the domain. 

In petroleum reservoir engineering and neutron diffusion problems quite often equations 
with strongly discontinuous coefficients appear. For these problems equations (4.5.3) and 
(4.5.4) are not representative. Suitable test problems with strongly discontinuous coefficients 
a\e been proposed in [111] and [79]; a definition of these test problems may also be found 
in [80]. In Kershaw’s problem the domain is non-rectangular, but is a rectangular polygon. 
The matrix for both problems is symmetric positive definite. With vertex-centered multigrid, 
operator-dependent transfer operators have to be used, of course. 

The four test problems just mentioned, i.e. (4.5.3), (4.5.4) and the problems of Stone 
and Kershaw, are gaining acceptance among conjugate gradient and multigrid practitioners 
as standard test problems. Given these test problems, the dilemma of robustness versus 
efficiency presents itself. Should one try to devise a single code to handle all problems (ro- 
bustness), or develop codes that handle only a subset, but do so more efficiently than a robust 
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code 7 This dilemma is not novel, and just as in other parts of numerical mathematics, we 
expect that both approaches will be fruitful, and no single ’best’ code will emerge. 


Numerical experiments for the test problems of Stone and Kershaw and equations 4.5 3) 
and (4.5.4), comparing CGS and multigrid, are described in [107], using LU and TOLU pre- 
conditioning and smoothing. As expected, the rate of convergence of multigrid is unaffected 
when the mesh size is decreased, whereas CGS slow down. On a 65 x 65 grid there ,s not great 
difference in efficiency. Another comparison of conjugate gradients and multigrid is prese 
in [40]. Robustness and efficiency of conjugate gradient and multigrid methods are deter- 
mined to a large extent by the preconditioning and the smoothing method respectively. The 
smoothing methods that were found to be robust on the basis of Fourier smoothing analysis in 
Chapter 4 suffice, also as preconditioners. It may be concluded that for medium-sized linear 
problems conjugate gradient methods are about equally efficient as mult, grid in acce teratmg 
basic iterative methods. As such they are limited to linear problems, unlike multigrid. On 
the other hand, conjugate gradient methods are much easier to program, especially when the 
computational grid is non-rectangular. 


7 Finite volume discretization 

7.1 Introduction 

In this chapter some essentials of finite volume discretization of partial differential equations 
are summarised. For a more complete elementary introduction, see for example [46] or [89]. 
We wiU pay special attention to the handling of discontinuous coefficients, because there seem 
to be no texts giving a comprehensive account of discretization methods for this situation. 
Discontinuous coefficients arise in important application areas, such as porous media flows 
(reservoir engineering), and require special treatment in the multigrid context. Furthermore, 
hyperbolic systems will be briefly discussed. 

7.2 An elliptic equation 

Cartesian tensor notation is used with conventional summation over repeated Greek subscripts 
(not over Latin subscripts). Greek subscripts stand for dimension indices and have range 1, 2, 
..., d with d the number of space dimensions. The subscript denotes the partial derivative 

with respect to x a . 

The general single second-order elliptic equation can be written as 


Lu = —(a a 0 u ,a),0 + (bail), a + cu = s in Q, C R d 


( 7 . 2 . 1 ) 
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The diffusion tensor a Q0 is assumed to be symmetric: a a0 = a 0a . The boundary conditions 

will be discussed later. Uniform ellipticity is assumed: there exists a constant C > 0 such 
that 

a a pv a vp > Cv a v a , Vt? G lR d (7.2.2) 

For d— 2 this is equivalent to Equation (7.2.9). 

The domain ft 

The domain ft is taken to be the d-dimensional unit cube. This greatly simplifies the con- 
struction of the various grids and the transfer operators between them, used in multigrid 
In practice, multigrid for finite difference and finite volume discretization can in principe be 
applied to more general domains, but the description of the method becomes complicated, 
and general domains will not be discussed here. This is not a serious limitation, because the 
current main trend in grid generation consists of decomposition of the physical domain in 
subdomains, each of which is mapped onto a cubic computational domain. In general, such 
mappings change the coefficients in (7.2.1). As a result, special properties, such as separa- 
bility or the coefficients being constant, may be lost, but this does not seriously hamper the 
application of multigrid, because this approach is applicable to (7.2.1) in its genera] form. 
This is one of the strengths of multigrid as compared with older methods. 

The weak formulation 

Assume that a is discontinuous along some manifold T C ft, which we will call an interface; 
then Equation (7.2.1) now has to be interpreted in the weak sense , as follows. From (7.2.1) 
it follows that 

(Lu,v) = (s,v) Vv£H, (u,v) = J uvdil (7.2.3) 

n 

where H is a suitable Sobolev space. Define 

a(u,v) = f a a0 u >oi v t0 dfl — f a a0 u a n 0 vdY 

Cl an ’ 

b(u,v) = f(b a u) <a vdfi. (7.2.4) 

n 

with n 0 the x 0 component of the outward unit normal on the boundary 0ft of ft. Application 
of the Gauss divergence theorem gives 

(Lu, v) = a(u, v) + b(u, v ) + ( cu , v) 

The weak formulation of (7.2.1) is 

Find ueH such that a(«,u) + b(u,v) + (cu,v) = (*,»), Vv <E H 


(7.2.5) 

(7.2.6) 
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For suitable choices of H, H and boundary conditions, existence and uniqueness of the solution 
of (7.2.6) has been established. For more details on the weak formulation (not needed here), 
see for example [31] or [58]. 


The jump condition , r 

Consider the case with one interface I\ which divides ft in two parts and Jl 2 , m each ot 
which a a n is continuous. At T,a a 0 (x) is discontinuous. Let indices 1 and 2 denote quantities 
on T at the side of fii and fi 2 , respectively. Application of the Gauss divergence theorem to 
(7.2.5) gives, if u is smooth enough in fl 1 and ff , 

a(u,v)=- J {a a 0 u ,a), 0 v ^ + J ( a a 0 U ]a ~ a a 0 V \a) n 0 V< ^ (7.2.7) 

o\r r 


Hence, the solution of (7.2.6), if it is smooth enough in Sli and ft 2 , satisfies (7.2.1) m fi \ I\ 
together with the following jump condition on the interface T 


a aP U ]a n 0 ~ a a0 V \a Tl P 


on 


(7.2.8) 


This means that where a a0 is discontinuous, so is u <a . This has to be taken into account in 
constructing discrete approximations. 

Exercise 3.2.1. Show that in two dimensions Equation (7.2.2) is equivalent to 


®ll a 22 — a 12 0 


(7.2.9) 


7.3 A one-dimensional example 

The basic ideas of finite difference and finite volume discretization taking discontinuities in 
a a p into account will be explained for the following example 

— (au^),! = 5 , x £ fl = (0, 1) (7.3.1) 

Boundary conditions will be given later. 

Finite difference discretization 

A computational grid G C ^ is defined by 

G = {x € M : x = Xj = jh, j = 0,1,2 ,...,n, h = 1/n} (7.3.2) 

Forward and backward difference operators are defined by 

A Uj = (tij+j - Uj)/h, Vuj - ( Uj - Uj-\)/h (7.3.3) 
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A finite difference approximation of (7.3.1) is 
symmetric formula is 


obtained by replacing d/dx by A or V. A nice 


2 (V(aA) + A(aV)}n^ = Sj , j = 1,2 1 


(7.3.4) 


where sj - a(xj) and Uj is the numerical approximation of u( Xj ). Written out in full 
Equation (7.3.4) gives ’ 


{-(a,., + a,K_i + ( aj _, + 2 a, + a j+1 ) Uj - (a, + a j+1 )u j+l }/2h 2 = Sj , 


j = 1,2 1 


(7.3.5) 


If the boundary condition at z = 0 is „(0) = / (Dirichlet), we elimine u 0 from (7.3.5) with 
u° - f If the boundary condition is «(0)«,,(0) = / (Neumann), we write down (7.3.5) 
for j = 0 and replace the quantity -(a., + + («_, + a 0 )u 0 by 2/. If the boundary 

condition is cittj(O) + c 2 u(0) = / (Robin), we again write down (7.3.5) for j = 0, and replace 

the quantity just mentioned by 2(f - c 2 u o )a(0)/c v The boundary condition at * = 1 is 
Handled in a similar way. 

An interface problem 

example that ^' 3 '^ ^ inaccurate for interface problems, we consider the following 

a(x) = f, 0 < x < x*, a(x)= 1, x* < x < l (7.3.6) 

The boundary conditions are: u(0) = 0, u(l) = 1. The jump condition (7.2.8) becomes 


e lim u j = lim u i 

x]x* x{x* 1 


(7.3.7) 


By postulating a piecewise linear solution the solution of (7.3.1) and (7.3.7) is found to be 

u = ax, 0 < x < x*, u = eax + l-ea, x* < x < 1 
a = l/(x* - ex* + e) ’ (7.3.8) 

Assume x* < x* < x k+i . By postulation a piecewise linear solution 

Uj-aj, 0 < j < k, uj = 0j-0n+ 1 , k + 1 < j < n ( 7 . 3 . 9 ) 

(7.3 9) trith th<! S ° IUti °" ° f <7 3 ' 5) ’ Witl1 ‘ he bo ’ , " dar> ' co “ di,i «" 5 e^en above, is given by 


Hence 


{3 = ea, a= ^ + e(n - k) + kj 


u k 


x k 


£ h( 1 - 0/(1 + s) + (1 — e)xk + e 


(7.3.10) 

(7.3.11) 
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(7.3.12) 


Let x* = Xk+\. The exact solution in x k is 

u{xk) = 

Hence, the error satisfies 


Xk 


(1 - + £ 


(7.3.13) 


u k - u(*fc) = o 

As another example, let *> = x„ + h/2. The numerical solution in x t is stiU given by (7.3.11). 
The exact solution in x k is 

x k (7.3.14) 


u ( x k) (\-s)xk + £ + h( l-e)/2 


The error in x k satisfies 


u k 


- = 0 (iiiTi)' 1 ) 


(7.3.15) 


When a(x) is continuous (e = 1) the error is zero. For general continuous a(x) the error is 
0(h 2 ). When a(z) is discontinuous, the error of (7.3.4) increases to 0{ ). 

Finite volume discretization . 

By starting from the weak formulation (7.2.6) and using finite volume discretization one may 

obtain O(f’) accuracy for discontinuous o(x). The domain fi is (almost) covered by cells or 
finite volumes Slj, ^ = ^ ^ + ^ = lj2 , , (7.3.16) 

Let w(z) be the characteristic function of 9,j 

v{x) = 0, x i «( x ) = !. x e (7.3.17) 

A convenient unified treatment of both cases: a(x) continuous and o(z) discontinuous is as 
follows We approximate a(x) by a piecewise constant function that has a constant value a 
in each Sl r Of course, this works best if discontinuities of »(*) he at boundaries of fin 
volumes flj. One may take dj — a(xj), or 


a 3 = h 1 j 


ad£t 


With this approximation of a(x) and v according to (7.3.17) one obtains from (7.2.7) 

a(u,u) = — J {au t \) t idil 

= — au,i fj+ h h % if 1 < J < « - 1 (7.3.18) 
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By taking successively j = 1 
unknowns Uj(u 0 = 0 and u n 


,2, n 1, Equation (7.2.6) leads to n - 1 equations for the n- 1 
— 1 are given), after making further approximations in (7.3.18). 


In ; rd i?rr oxima,e auMxi + hm we proc “ d •* ««. is smooth 

u'tx t S d r, n0 j°° !l W j " mPS 11 1 = X ’ +h/i - Hen “' M is a bad idea *« discretize 

u A\ x j + h/2), Instead, we write 


X 3 + l 


r j+ 1 


-f-i f t 1 

=7 u ’' dx = j -« w ,i^ = (a«i) J+ i /2 [ -dx 

x > *, 4 a 

where we have exploited the smoothness of au A . We have 

x i + i 

/ ^ = h ! w i 

x i 

with Wj the harmonic average of a, and Oj+j: 


tn, 


— 2ajaj + i/(aj + aj +1 ) 


(7.3.19) 


(7.3.20) 


(7.3.21) 


and we obtain the following approximation: 

( au ,i)j+i /2 - Wj(uj +1 - Uj)/h (7.3.22) 

^creUzaUon nS (? ' 3 ' 18) and (7 ' 3 ' 22) ’ the weak Emulation (7.2.6) leads to the foUowing 


- Wj-O/h - uij(ti j+1 - Uj )/h = hsj, J r= 1,2 , ..., n — 1 


(7.3.23) 


with 


--7 

«> 


sdz . 


(T.3“)“ (X) ^ Sm ° 0th ' W ’ ” ^ + °' +1,/2 ’ and We reC ° Ver the * nite differe "“ approximation 

Equation (7.3.23) can be solved in a similar way as (7.3.5) for the interface problem under 
consideration. Assume Z* = x k + h/2. Hence 


w i -e, 1 <3 <k; w k = 2e/(l + e); Wj = 1, * < j < „ _ i 

Again postulating a solution as in (7.3.9) one finds 

P = ae, a = w/[e - we(k + 1 - n) + wJk] 


(7.3.24) 

(7.3.25) 
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a _[(l_ e )/2 + e(n-k) + k]~' =h/[{x k + h/2){\-e) + e] (7.3.26) 

Comparison with (7.3.8) shows that u 3 = «(*;): the numerical error is zero. In more general 
circumstances the error will be 0(h 2 ). Hence, finite volume discretization is more accurate 
than finite difference discretization for interface problems. 

Exercise 3.3.1 The discrete maximum and I 2 norms are defined by, respectively, 

, 1/2 


Itt 


OO 


= maxljuj) : 0 < j < 


n}, 


(7.3.27) 


Estimate the error in the numerical solution given by (7.3.9) in these norms. 

7.4 Cell-centered discretization in two dimensions 

TheXmtin^ if divided in cells as before, but now the grid points are the centers of the 
cells, see Figure 7.4.1. The computational grid G is defined by 


1 1 . 

G = {x G ft : x = x 3 = (j - s)h, j = (ji,ja), * - V' 

h = (hiih?), j a = 1,2, h a — 1 ! n ot) 


(7.4.1) 


The cell with centre x 0 is caUed Qj. Note that in a ceU-centered grid there are no grid points 
on the boundary dSl. 

Finite volume discretization in two dimensions 

Integration of (7.2.1) over a finite volume gives, with c = s = 0 for brevity, 


— J a a pU' a n0dT + J b a un a dT — 0 


(7.4.2) 


with T. the boundary of ft, and n the outward unit normal. Let us (denote) the ’’east” part 
of Tj, at x\ = (ji + by IV On IV n = (1,0)- 


The convection term 
We write 


J b\udT = h 2 (b\u)j 1+ \/ 2 ,j 2 


(7.4.3) 
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Figure 7.4.1: Cell-centered grid. (• grid points; finite volume boundaries.) 


Central discretization gives 

( b l u )ji+l/2j, - \{{ b \u)j + (biu) jl+hh } (7.4.4) 

Upwind discretization gives 


( 6 i u )ii+i/2,i, - ^{(bi + |6i|ti>j + ^{(6j - |MMji+i,j 2 (7.4.5) 

If a 12 = 0 then (7.4.4) results in a A'-matrix (definition 3.2.5) only if with w defined below 


^I^ji+i»j2 1 2 I 2 

W h+i/2,h ~ ’ W ji-l/2,h ~ 


(7.4.6) 


whereas, if a 12 = 0, (7.4.5) always results in a A'-matrix. The advantages of having a A'-matrix 
are 


• Monotonicity: absence of numerical wiggles. 

• Good behaviour of iterative solution methods, including multigrid, as discussed in Chap- 
ter 4. 

The diffusion term 
We write 

J a al u tC ,dr “ h 2 (a al u^) ji+l/2 j2 (7.4.7) 

r e 
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and approximate the flux Fj i+ 1 / 2 , j 3 
dimensional case, taking into account 


= (a al u a ) n+1/2lj2 in a similar way as in the one- 
the fact that (u,\)j+\j 2 ,i 2 does not exist, but F and u )2 


are smooth. We have 


ji+i J2 


u\ n ^ 

U huj2 


ar Jl+ 1, J2 

f -^dx 1 = hi Wj 

J au 

Xj 


X h+l>J2 X J1+1.J2 

f U'ldx! = / 

x 3 1<32 


* 31 <32 


*21 + 1.1/2 

f ^r{( a ia u .«) - a i2 w ,2)aa:i 

n ’ n *>1+1/2, j 2 Xj i + 1 ' 1/2 , , 

S F jl+1/2ij2 / ^7^1 - (*, 2)j,+i/2,j2 / “n/audx. 


(7.4.8) 


We now assume that ai 2 = 0, or else that u )2 may be approximated by 


(n, 2 )j,+i/2,j2 - 4/^“ 


jl,>2+l , |il+l>J2+ll 


(7.4.9) 


This will be accurate only of a a p is smooth at the north and south edges, the general case 
seems not to have been investigated. 


We obtain, with 


X J1 + 1 .J2 

u, ii+i/2,j2 - fe / / 

X JlJ2 


(7.4.10) 


r Jl + 1 .J2 

With (7.4.9) the if-matrix property is lost. The off-diagonal elements with the wrong side 
are much smaller than those generated by central discretization of the convection term at 
high Peclet number, and usually results obtained and performance of iterative methods are 
still satisfactory. See [141] for more details, including a discretization that gives a A -matrix 

for ai 2 ^ 0. 


7.5 A hyperbolic system 
Hyperbolic system of conservation laws 

In this section we consider the following hyperbolic system of conservation laws : 

^ + + = (*,„)€ n, t e (o,t] (7.5.1) 

dt dx dy 
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where 


w:[0,T]xft^S a cM, s :[0,T]x ft - 1W > , f,g-.S a ^BV> (7.5.2) 

Here S a is the set of admissible states. For example, if one of the p unknowns, u, say, is 
the fluid density or the speed of sound in a fluid mechanics application, then u, < 0 is not 
admissible. Equation (7.5.1) is a system of p equations with p unknowns. Here we abandon 
Cartesian tensor notation for the more convenient notation above. Equation (7.5.1) is assumed 
to be hyperbolic . 

Definition 7.5.1 Equation (7.5.1) is called hyperbolic with respect to t if there exist for all 

¥> € [0, 2tt) and admissible u a real diagonal matrix D(u,(p) and non-singular matrix R(u <p) 
such that ’ 

A(u, ip)R{u, <p) = R(u , <p)D(u, (p) (7.5.3) 

where 

a( \ df( u ) . dg(u) 

A(u, (p) = cos + sin y>- £ ; (7.5.4) 

The main example to date of systems of type (7.5.1) to which multigrid methods have been 
applied successfully are the Euler equations of gas dynamics. See [34] for more details on 
the mathematical properties of these equations and of hyperbolic systems in general. For 
numerical aspects of hyperbolic systems, see [101] or [104]. 

For the discretization of (7.5.1), schemes of Law-WendrofF type (see [101]) have long been pop- 
ular and still are widely used. These schemes are explicit and, for time-dependent problems, 
there is no need for multigrid: stability and accuracy restrictions on the time step At are 
about equaly severe. If the time-dependent formulation is used solely as a means to compute 
a steady state, then one would like to be unrestricted in the choice of At ands/or use artificial 
means to get rid of the transients quickly. 

In [92] a method has been proposed to do this using multiple grids. This method has 
been developed further in [76], [30] and [77], The method is restricted to Lax-Wendroff type 
formulations. 

Finite volume discretization 

Following the main trend in contemporary computational fluid dynamics, we discuss only the 
cell-centered case. The grid is given in Figure 7.4.1. Integration of (7.5.1) over gives, using 
the Gauss divergence theorem, 

— J udQ + J ( f(u)n x + g{u)n y )dY = f Sdtl (7.5.5) 

r j i 
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(7.5.6) 


where is the boundary of Sl } . With the approximations 

j udil ~ \Qj\uj, J SdSl ~ \flj\sj 
u, n, 


where |fij| is the area of ftj, Equation (7.5.5) becomes 

\Slj\duj/ dt + J ( f{u)n x + g(tt)n v )dr = 


(7.5.7) 


The time discretization will not be discussed. The space discretization takes place by approx- 
imating the integral over r ; . 

Let A = Xj + {h\l% -fca/2), B = xj + (hi/2,h 2 /2), so that AB is part of E,. On AB, n x = 1 
and n v = 0. We write 

B 

J f{u)dx 2 3 h 2 /(tt) c 

with C the midpoint of T5. Central space discretization is obtained with 


(7.5.8) 


f(u)c = \ f(Uj) + \f( U J+e l) 


(7.5.9) 


In the presence of shocks, this does not lead to the correct weak solution, unless thermo- 
dynamic irreversibility is enforced. This may be done by introducing artificial viscosity, an 
approach followed in [72]. Another approach is to use upwind space discretization, obtained 

by flux splitting: 

f(u) = f + (u ) + / (u) (7.5.10) 

with /*(«) choosen such that the eigenvalues of the Jacobians of /*(«) satisfy 


A (df + /du) > 0, A (df~/du) < 0 


(7.5.11) 


There are many splittings satisfying (7.5.11). For a survey of flux splitting, see [64] and [126]. 
With upwind discretization, f(u)c is approximated by 


f(u)n a f + (u, + f-ftti+e, )) 


(7 R 19.1 


The implementation of boundary conditions for hyperbolic systems is not simple, and will 
not be discussed here; the reader is referred to the literature mentioned above. 

Exercise 7.5.1 Show that the flux splitting (7.4.5) satisfies (7.5.11). 
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8 Conclusion 


An introduction has been presented to the application of multigrid methods to the numerical 
solution of elliptic and hyperbolic partial differential equations. 

Because robustness is stongly influenced by the smoothing method used, much attention has 
been given to smoothing analysis, and many possible smoothing methods have been presented. 

An attempt has been made to review much of the literature, to help the reader to find his 
way quickly to material relevant to his interests. For more information, see [141]. 

In this book application of multigrid to the eqations of fluid dynamics is reviewed, a topic 
not covered here. There the full potential equation, the Euler equations, the compressible 
Navier-Stokes equations and the incompressible Navier-Stokes and Boussinesq equations are 
treated. 

The principles discussed in these notes hold quite generally, making solution possible at a cost 
of a few work units, as discussed in chapter 6, for problems more difficult than considered 
here. 
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